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CALCULATE THE STRESSES OF ONE CRACK 
EMBEDDED IN ISOTROPIC PLATE 

This program coded in Macsyma calculates the stresses of a crack embedded 
in isotropic plate by solving the governing equation with certain boundary conditions. 

1. Define the governing equation 

d A 2d 3 4 8 4 

^ F(x ' y)+ W^y F(x ' y)+ ¥-y F{x ’ y) = l1 ’ (1) 

where F(x , y ) is the stress function. 

By Fourier transform the governing equation can be changed to an ODE 
with constant coefficients respect to y. 

If F(x,y ) is evaluated by J^° 00 (f>(s,y)exp~ t,x ds the integration will be 
performed by Macsyma. Since only the integrand is involved during the calculation, 
we can omit the integral sign. The integral sign will be added later on. 

2. Taking the Fourier transform yields the ODE 

s A 4> ~ s 7 4>" + 4>"" = 0 . ( 2 ) 

3. Evaluating <f> by C exp ry * in equation (2) the characteristic equation is obtained 

r 4 — 2r 2 + 1 = 0. (3) 


1 
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4. Solve the characteristic equation. 

5. If the roots of the characteristic equation are two identical real roots, the general 
solution of equation (2) is 

(C[ + C' 2 y) exp 1 " +(C' 3 + C[y) exp” 1 " . (4) 


6. Since the stress function F(x,y) is bounded at infinity, the solution of the gov- 
erning equation (1) is 


F(x,y + ) = + C 2 y)exp |,|5/ exp ” x ], 

for y > 0, 

(5) 

F (*, y ~) = ^[(<?3 + C 4 y) exp |j|v exp -,,x ], 

for y < 0, 

(6) 


where C,-(t = 1, ...,4) are arbitrary functions of s. 

Applying the continuity conditions Cz and C 4 can be expressed in terms of 
G\ and C 2 . 


7. Define equation 


^(z,0+)-F(z,0-) = 0 


and by evaluating F(x,y + ) and F(x,y ) at 0,we obtain 


Cz — C\. 


8. Define equation 

dF(x, 0 + ) dF(x,0~) _ Q 

dy dy 


( 7 ) 

(8) 


( 9 ) 
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by differentiating equuation (9) w.r.t. y and evaluating y at 0. Solving C\ and get 


C 4 = -2|s|6i + c 2 . 


Displacements u and v are obtained from Hooke’s law 


( 10 ) 


v — d - an J F xi ;dy i 

u = a n J Fyydx + a n F x , 


( 11 ) 

( 12 ) 


where an and a 4 2 are coefficients of plane stress and plane strain. 

9. Computing u + and u~ by Macsyma. 

10. Computing v + and v~ by Macsyma. 

If we introduce two auxiliary functions fi(x) and f 2 (x) properly, C 2 and C 2 
can be written in terms of fi(x) and f 2 (x). 

11. Define 

fi( x ) = ^[« + (*i°) •“«"(*» °)1- (13) 

After differentiating u + and u~ w.r.t. x and evaluating y at 0, f\(x) is simplified by 
Macsyma simplication functions. So 

(2a.iiC72 1^| — aanCiS 2 ) exp -WI 


M x ) 


7T 


(14) 


12. Define 


/*(*) = ^[ v+ ( z > 0 ) “ v ( z »°)]- (15) 

After differentiating v + and v~ w.r.t. x and evaluating y at 0, / 2 (z) can not be 
simplified only by Macsyma simplification functions. There are terms such as 


g(x, s)s 2 + h(x , s)s 3 -f tu(x,.s).s|s| 


(16) 


and these terms can be simplified by substitutions. For example, s 3 is substituted 
by &|.s|, then k is substituted by j|j|. Although it is the same as .s 3 = ,s|,s||.s|, the 
only difference is that terms like -jyj are simplified to s|s|. Finally we obtain 


M x ) 


2ia 11 C' 1 s|s| exp - ” 1 
7 r 


(17) 


13. Taking the inverse Fourier transform on both sides of equation (15) (Notice that 
we omitted the integral sign during the calculation). C\ is expressed in terms of / 2 (x), 

C\ = V, fz{t) exp”‘ . (18) 

4a n js|s 

14. Taking the inverse Fourier transform on both sides of equation (14) (Notice that 
we omitted the integral sign during the calculation). With equation (18) C 2 is ex- 
pressed in terms of /i(z) and / 2 (z), 
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15. Calculating C 3 from equation ( 8 ), we have 


C 2 — 


4a.ii 


/ 2 (t)exp w . 


( 20 ) 


16. Calculating C 4 from equation ( 10 ), we have 


c 4 = -i-[i/,(f)exp i “+-/ ! (()exp i "]. 
4a n \s\ 3 

( 21 ) 

17. Stresses cr xx cr y y and r xy are defined as follows 

d 2 F(z,y) 

“ a-y ’ 

( 22 ) 

“ 5 2 x ’ 

(23) 


(24) 


18. Simplification of a xx . 

A second derivative of F(x,y) w.r.t. y and substitution of Ci(i = 1 , ..., 4) 

yields 

is f2(t)y exp~’ v+ ' ,x ~' ,t _ sfi(t)y exp~ ,y+ ' ,x ~ ,,t 
4 a n 4 a n 

, if 2 (t)exp- ,v+i ' x -' ,t , fjjt) exp-' v+ijs ~ ijt 
4an 2an 

isf 2 (t)y exp -'v-wx+i.t sfi(t)y exp~‘ y ~' ,x+iat 

4au 4au 
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i/ 2 (t)ex p -«v-'«+»* /i(<)exp- ,y -‘”+» rt 


4a 


11 


2a 


li 


(25) 


With integral signs, we know the stress a sx is 



isf 2 (t)y exp~* y+ ** z ~ 1 * t 
4a n 


sfi(t)yexp ’v+ i,x l,t 

4 a u 


| i'/ 2 (f ) exp~* v+t * J ~** t /i(f)exp-» +, '"- , ' ,t 
4a 2 j 2au 

isf 2 (t)y exp~ ,v_wz+t,t s f\{t)y exp~* v ~” :E+, * t 
4an 4an 


t7 a (t)exp-*v- f “^ [ ACQexp-*^”*"* 


4a 


li 


2a 


ii 


(26) 


The function integrate in Macsyma can not be used directly to to perform the 
integration of equation (25) w.r.t. 5, because it keeps calculating and never comes 
to an end. Also we know the fact 


f exp^ a±w )* ds = 
Jo 

1 

(27) 

a ± bi ’ 

[ sexp^ a±w ^<fs = 

JO 

1 

(28) 

(a ± bi) 2 


Substitution °f — f° r all expt v±, ^ z terms and p^X(~ 7 )jp for all y exp^ v ~^ x £ )b 

terms has the same effect as integration. So the following expression is obtained 

1 r KfzCQ ~ h{t))y f ifijt) hit) 

87ra!i [*(f - x) - y} 2 [~y - i{t - x)] 2 [-y - i(t - z)p J 
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2/i(t) - i/ 2 (Q if 2 {t) + 2fi{t) 
i(t - x) - y -y - i(i - x) 


}• 


( 29 ) 


Taking the common denominator with the same order of power, we ultimately obtain 


= 


47ra 


u 


-{/*(*)(*-*)[ 


-2 y 2 


+ 


- y/i(0[ 


[y 2 + (i - x) 2 ] 2 y 2 + (t — x) 2 
y 2 -{t-x) 2 2 


[y 2 + {t - x) 2 ] 2 y 2 + (t - x) 2 


]}• 


(30) 


Finally by adding the integral sign, we have the exact answer 


- ik }> m - x)[ VW^n + 

y 2 - (t - x) 2 2 


-y/i(0[ 


[y 2 + (£ - x) 2 ] 2 y 2 + (t- x) 2 


}}dt. 


(31) 


19. Following the same procedures to simplify cr w with slight change, we obtain 


'w 


= 4^1 - l)[ 


2y 2 


+ 


+ y/i(0 


[y 2 + (t-x) 2 ] 2 y 2 + (t — x) 2 
y 2 -(t-x ) 2 


) 


[y 2 + (t - x ) 2 ] 2 


}dt. 


(32) 


20. Following the same procedures to simplify r xy with slight change, we obtain 


-1 r 


Tiy “ 4™ 


y 2 -(t- x) 2 


y 2 + (t- x) 2 


- vh{t) 


[y 2 + (t - x) 2 } 2 


}dt. 


(33) 



CALCULATE THE STRESSES OF TWO CRACKS 
EMBEDDED IN ISOTROPIC PLATE 

According to local coordinate system the stresses for one crack is already 
given. For two cracks, the relation between the local coordinate systems x x ,y x and 
x 2 ,y 2 are 

r lx + X\ cos ifi - y x sin <p x = r 2x -f x 2 cos y> 2 — y 2 sin (p 2 (1) 

r ly + x x sin <p x + y x cos = r 2y -f x 2 sin <p 2 + y 2 cos <p 2 (2) 

the transformation equations for plane stress are 

cr X2XJ = ex xjxi COS 2 9 + a yiVl sin 2 9 + r XlVl sin 29 (3) 

a v2V2 - s ^ 2 0 + a vivi cos2 & ~ T *ivi s5n ( 4 ) 

r x 2 V 2 = — °vm] s * n 20 + t XiVi cos 29 (5) 

where 9 = (p 2 — y x . 

Thus, the total stress for the first crack is the summation of the stress of 
the first crack and the transformed stress from the second crack. The total stress 
of the second crack is vice versa. Let x,- = a,£ and U = a,r , where i = (1,2), by 
applying the boundary conditions, the total stress for crack 1 can be expressed 

< = £ 7=i iT + £ k * fniT + £ k iV» iT y < 6 ) 
^ £ r~i dT + £ k » h ' iT + £ k %f” iT } ( ? > 


i 
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where k are Fredholm kernels (a,/3 y n,m = 1,2) 


jj 0,2 , sin 2$(p 4 — a 3 £sin 0)[(a 2 r — p 3 — a 3 £ cos 0) 2 — ( p 4 — a 3 £sin 0) 2 ] 

21 7T [(p 4 — a^sin 0) 2 + (a 2 r — p 3 — ai^cos ^) 2 ] 2 

2cos 20(a 2 r — p 3 — a 3 £cos fl)(p 4 — a 3 £sin 0) 2 
[(p 4 - a^sin Oj 2 -f (<z 2 r - p 3 - Gi<fcos 0) 2 ] 2 


^ sin 2$(p 4 — a 3 £cos 6 ) — cos 2#(a 2 r — p 3 — a^cos 0) 
(p 4 — a 3 £sin 0) 2 + (a 2 r — p 3 — a 3 £cos0) 2 ' 


( 8 ) 


. u _ a 2 r cos 20(p 4 — ai^sin 0)[(a 2 r — p 3 — a^cos 0) 2 — (p 4 — aif sin 0) 2 ] 
22 7T [(p 4 — a^sin 9 ) 2 + (a 2 r — p 3 — a a £cos 0) 2 ] 2 


2sin 20(a 2 r — p 3 — a 3 £cos 0)(p 4 — ai£sin 6 ) 2 
[(p 4 — a 2 £sin 0) 2 + (a 2 r - p 3 — a^cos 0) 2 ] 2 


( 9 ) 


^.12 _ a 2 r 2sin 29(a 2 r — p 3 — a^cos 0)(p 4 — ai£sin 0) 2 
21 ^ [(p 4 — ai£sin 0) 2 -f (a 2 r — p 3 — a 3 £cos 0) 2 ] 2 

cos 2#(p 4 — a^sin 9)[{a 2 T — p 3 — aj£cos 0) 2 — (p 4 — aj.£sin #) 2 ] 
[(?4 ~ a^sin 0) 2 + (a 2 r - p 3 - cos 0) 2 ] 2 


^ 2sin 2 0(p 4 — a^cos 0) — sin 20(a 2 r — p 3 — ai<fcos 0)^ 
(p 4 — a^sin 6 ) 2 + (a 2 r — p 3 — ai£cos0) 2 


( 10 ) 
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£12 _ a 2 rSin20(p 4 — a x £sin 0)[(a 2 r — Pz — a x ( cos 0) 2 - (p 4 — a x £sin 0) 2 ] 
7T [(p 4 — a x £sin 0) 2 + (a 2 r — p 3 — ai£cos #) 2 ] 2 


2cos 2#(a 2 r — p 3 — a x £cos #)(p 4 — a^sin 6) 2 
[(p 4 - aj^sin 0) 2 + (a 2 r - p 3 - a x £cos 0) 2 ] 2 


where 


+ 


(a 2 r -p 3 - a^cosO ) 


(p 4 — ai^sin 0) + (a 2 r — p 3 — a x £co.s0) 2 


} 


( 11 ) 


£3 = (ri y - r 2 „) sin p 2 + (r lx - r 2x ) cos ip 2 


Vi = ( r iv - r 2y) cos <p 2 - (r lx - r 2x ) sin <p 2 


The total stress for crack 2 can be expressed 

< = i^ { ; L rt\ iT + L k " f ' idT + L kllh ' iT} (12) 

< t2 w = 7—{-[ ——zdT + f kllf u dT+f k 2 lf i2 dr} (13) 

4a lx 7r J- i r — £ j-i J-i 


where are Fredholm kernels (a,/3,n,m = 1,2) 


£21 = £l{_^ 


2 0(p 2 + a 2 <fsin #)[(a x r — p x — a 2 <fcos #) 2 — (p 2 + a 2 £sin #) 2 ] 
[(p 2 + a 2 f sin 0) 2 + (air - p x - a 2 £cos 0) 2 ] 2 


2cos 2#(a x r — p x — a 2 £cos #)(p 2 + a 2 £sin 6 ) 2 
[(p 2 + a 2 £sin 0) 2 + (air - p x — a 2 fcos 0) 2 ] 2 


+ 


— sin 


20(p 2 + a 2 £sin 6) — cos 2#(a x r — p x — a 2 £cos #) 
(p 2 + a 2 £sin 0) 2 + (a x r — p x — a 2 (cos0) 2 


} 


( 14 ) 
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, 21 _ a i / cos 2#(p 2 + a 2 £sin 0)[(air — pi — a 2 ( cos #) 2 - (p 2 + a 2 £sin 0) 2 ] 
12 [(p 2 + a 2 ^sin 0) 2 + (air — pi — a 2 ^cos ^) 2 ] 2 

2sin 20(ajr — pi — a 2 £cos 0)(p 2 + a 2 £sin G) 2 
[( P 2 + a 2 £sin G) 2 -f (air - p : - a 2 (cos G) 2 ] 2 


(15) 


^ 22 _ a i / — 2sin 20(air — pi — a 2 £cos 0)(p 2 + a 2 £sin 0) 2 
U T [(?2 + ^sinfl) 2 + (air — p x — a 2 £cos 0) 2 ] 2 

cos 2$(p 2 + a 2 (sin 0)[(air — pi — a 2 fcos 0) 2 — (p 2 + a 2 fsin 0) 2 ] 
[(p 2 + a 2 £sin G) 2 + (air - Pi - a 2 ( cos 0) 2 ] 2 


2sin 2 0(p 2 + a 2 £sin 0) -f sin 20(a a r — pi — a 2 f cos #) 
(p 2 -f a 2 £sin G) 2 + (a 2 T — p\— a^cosO) 2 J 


(16) 


22 a x — sin 20(p 2 + a 2 ( sin 0)[(air — pi — a 2 £cos G) 2 — (p 2 + a 2 f sin 6 ) 2 ] 


k\l = -i{ 


[(p 2 + a 2 <fsin 0) + (a a r - pi - a 2 ( cos G) ] 2 


2cos 20(air — pi — a 2 (cos 0)(p 2 + a 2 £sin 0) 2 
[(p 2 + a 2 (sin G) 2 + (air - pi - a 2 £cos 0) 2 ] 2 


+ 


(air — pi — a 2 (cos 6) 


(p 2 -h a 2 £s in 6) + (air — pi — a 2 £cosG ) 2 


} 


( 17 ) 



where 
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Pi = ( r 2v _ T lv) Sin Pi + ( r 2* “ r l*) COS Pi 
Vi = ( r 2v “ r lv) C0S Pi — ( r 2* ~ r l*) sin Pi 



CALCULATION OF STRESS INTENSITY FACTORS 
OF n CRACKS EMBEDDED IN ISOTROPIC PLATE 



( 7 ) 
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r *V = ~ a w\ sin2(y?' - <p) + r xy cos 2(<p' - <p) (8) 

Label t' and a' are the stresses transformed by r < * and <x < < from x ,v' 

xy yy x y y y i & 

coordinate system to x,y coordinate system. Let x = a£, t = ar , and y = 0, where 
— 1 < (r < 1. After simplification we get 


K v = i k i(T^,r x ,r x ',r v ,r v >,<p,p,a,a')f' 1 dr 


+ / r l} ry , r y , ry , y?, y?', a, a')/'<fr} 

( 9 ) 

,a,a')f 1 dr 


r l , , 

+ J M r i£>W«SWyS^P> a » fl ')/i < * r } 

(10) 

Where 


Pi = (r y - r y < ) sin y>' + (r r - r s < ) cos y/ 

(11) 

P 2 = ~ V ) cos y?' - (r x - r x > ) sin y?' 

(12) 


and k; are Fredholm kernels ( i = 1, ...,4) and 6 = p — ip. 
ki(r,Z,r x ,r x ' , r„, r y < , y?, y?', a, a ') 


a . — sin 2 0(p 2 + a£sin 0)[(aV — p\ — a ( cos 8 ) 2 — ( p 2 + a£sin B) 2 } 
7T [(p 2 -f a£sin 0) 2 + (o'r — p\ — a£cos B) 2 } 2 

2c os 2B(a'r — pi — a£cos 8)(p 2 + afsin B) 2 
[(p 2 + as s i n ^) 2 + { a ' T — Pi — a(cos 8) 2 } 2 


— sin 2 #(p 2 4- a£ sin 8 ) — cos 2 #(aV — pi — a£cos #) 
(p 2 + a£sin 0) 2 + (aV — pi — a£cos8) 2 


( 13 ) 
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k 2 (r, t, r x , r x ., r„, r y . , <p, tp\ a, a') 


a f cos 26(p2 + a£sin 6 ){[ar — p\ — a£cos 6 ) 2 — (p2 + a£sin #) 2 ] 
7 T [(p 2 + a£sin 6 ) 2 -f (a'r — pi — a£ cos 0) 2 ] 2 


2sin 20(aV — pi — a£ cos 0)(p 2 + a£sin 0) 2 
[(p 2 + a£sin 0) 2 -(- (a'r — p\ — a £ cos #) 2 ] 2 


(14) 


*3(r, t,r x ,r x >,T y ,r v >,(p,<p',a, a) 


a . — 2sin 26(ar — p x — a£ cos 0)(p 2 + a£sin 0) 2 
7T [(p 2 4- a£sin #) 2 -f (a'r — pi — a£cos 0) 2 ] 2 

cos 20(p 2 + a£sin 0)[(a’r — Pi — a£cos 0) 2 — (p 2 + a£sin 6 ) 2 ] 
[(p 2 + a£sin 0) 2 4* (a'r ~ pi — a£cos 0) 2 ] 2 


+ 


2sin 2 0(p 2 + a£sin 0) -f sin 2 0(aV — pi — a£cos #) 
(j> 2 + a£sin #) 2 + (a'r — pi — a^cosB) 2 


} 


(15) 


ki( r > £, r„ r s . , r v , r » , y>, y?', a, a') 


a* — sin 20(p 2 + a£sin 0)[(aY — pi — a£cos 6 ) 2 — (p 2 + a£sin 0) 2 ] 
[(?2 + a<fsin 0) 2 + (a'r — p : — a£cos 0) 2 ] 2 

2cos 2 0(a'r — pi — a£cos #)(p 2 + afsin#) 2 
[(p 2 + a£sin 0) 2 -f (a'r — pi — a£cos #) 2 ] 2 

(a’r — p : — a£cos #) 

(p 2 + a£sin 6) 2 + (a'r — pi — a£cos6) 2 


( 16 ) 
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For two cracks, the total stresses for crack one are 


1 1 f 1 /n f 1 

' 1T ~ {“ / -dr + / k 1 (T,(,r Xl ,r X3 ,r Vl ,r V7 ,<p l ,<p 2 ,a 1 ,a 2 )f 21 dT 

7T J-l T — q V-l 

+ Li k ^ T} Tsi ’ 7-13 » r w » ’ V'l* > fll > a 2)/22^r} (17) 


T * v 4a 


u 


,ir_ 1 {I f JlL-dr+l k i (T,t,T Xl ,T X 7 ,r yn ,T Vi ,<p- L ,ip 2 ,a- L ,a 2 )f 2l dT 
ft J- i r — £ J- 1 

+ /_1 fc 4( T >^. r *,, r * J » r V l » r W,V :, 1 >V :, 2» a l. a 2)/22^'r} (18) 


<7 = 

w 4a 


li 


The total stresses for crack two are 


2T 

xy 


= 7— {- f “^7 dr+ [ *i(r,^,r^,r Xl ,r w ,r yi ,^ 2 ,y>i,o 2 ,ai)/ii<£r 
4an ft J -i T — q J - 1 

+ J fc 2 (t, £, r x , , r Xl , r w , r Vl , y> 2 , V>i , “2, ai)/i 2 <fr} (19) 


_ 2r _ {I / [ fcs( r »fi r xai r «ii r i8» r vx»V>2,¥»i,a2,oi)/ii^ 

ft J - i t — q J - 1 

r x 3 J r *: > r vj > r vi yf2,<Pi, dz, 0 ’\)f\zdr } (20) 


^ 4a 


li 


Where and r‘^ (z = 1, 2) should be equal to the fax field stresses trans- 
formed to each local coordinate system in opposite direction, 


C7 iT = 
w 


- & xx sin 2 ipi + a ly cos 2 ip { - r° v sin 2^,-j 


r iT = 
xy 


- 7') sin 2 f. + r ”» cos Ml 


(21) 

( 22 ) 
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So we get four integral equations which can be solved numerically. 

1 /-i f u /-i 

- 7 dr + h{T,t,r Xl ,r XJ ,r t/l ,r V 2 ,<p 1 ,<p 2 ,a 1 ,a 2 )f 2 idT 

7T J- 1 T — { J- 1 

+ k 2 (T,(,r Xl ,r X3 ,r Vl ,r V3 ,<p 1 ,<p2,a 1 ,a 2 )f22dT = 4a n r£ ( 23 ) 


1 /-l f 12 rl 

- 7 dr + fc3(T,^,r Xl ,r XJ ,r Vl ,r w ,V5 1 ,^ 2 ,ai,a 2 )/2icir 

7 r J- 1 t — £ y-i 

+ J_ i k 4 (T^,r Xl ,r Xi ,r yi ,r V3 ,ip 1 ,ip 2 ,a li a 2 )f22dT - 4a n cr^ ( 24 ) 


1 /■! / 21 /■! 

- 7 dr + ki(T,t,r mt ,r xli r ya ,T yit <p2 t <pi,a2,a l )f l idT 

7 r j-i t — 4 •'-l 

+ />•<• r i3> r sn r V3? r vn¥ > 2j¥ 5 ij02> a i)/i2<f r — 4aiiT xl/ (25) 


- / — - 2 - 2 — dr + / ^(t,^, r XJ j r*! ) r V3 J Tyx > ^2, Y?l, a 2j 0 -\)fl\dr 

7 r v-i t — £ J - 1 

+ /_ x fc 4(r,^r X3 ,r Xl ,r Vj ,r yi ,(p 2 ,^i,a2,a 1 )/ 1 2dr = 4a n cr^ (26) 


For n cracks, we get 2n integral equations 

-/ 7 dT + k \{r,i,r Xl ,T X7 ,---)f 2l dT + / k 2 (r,£,r Xl ,r XJ , • • -)f 22 dr 

7 r j-i r — J-i /-i 

• • + j ^ ki(r,{,r Xl ,r Xnf - • -)/m^ + * 2 (r,£,r Xl , r x „, • • -)/n 2 <fc- = 4a u r ;z 1 J (27) 

1 f 1 /12 7 1 7 1 

~ 7 dr + h(r,Z,r Xl ,r Xj ,---)f 2l dT + k 4 (r, (, r Xl , r Xj , • • -)f 22 dT 

7T y-1 T — £ y~i ^-1 

• + J ^h{r,Z,r Xl ,r Xn ,- ■ -)f nl dr + k 4 (r,(, r xx ,r Xn , • • -)/n2^ = 4a u cr^ (28) 
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J_ i h (r, r Ia ,r Xl , • • • )f u dr + J ^ k 2 (r,^r Xi , r Xl , • • -)/i 2 <*r + ^ ~r^ dr 
’• + J\ h(T,t,r X 3 ,r Xn ,-- .)UdT + k 2 (T,£,r X] ,r Xn ,- • -)/n 2 ^r = 4a u r x 2 J (29) 

J i k 3 (r,t,r X] ,r xll -..)fndT + k 4 (r,t,r X 3 ,r Xl , • • -)/i2^r + - yZT^ dr 

■■+ J i fc 3 ( r > £, r Xj , r x „, • • .)/-*■ + k<(T, (, r XJ , r Xn , • • -)/n 2 <£r = 40!!^ (30) 


’"in J T Xi 1 ‘ ‘ •)fn d T + J ^ k 2 {r,t, 

'/>< T 1 £ ) r *» > r *3 > ■ * •)f2\dr + J ^k 2 (r,(, 


T X n > r il 7 ' 


•)/l2^T 

. -)f 22 dr 


- f ~^~dr = 4a n r. 
7 r J-iT — t 


nT 

=v 


jHT,t,r Xn ,r Xl ,. • )fudr + J ^ k 4 (r,(,r Xn ,r Xli • • •)/«*- 


+ 


/.j * 3 (r,£,r Xn ,r X3l . • •)/«*■ + r x„ > r *2 ) ■ * ■)f22 d T 

... Lf'Js 

( 


- / ——7 dr = 4a n <7 ] 
ir J-i t — t 


nT 

w 


(31) 


(32) 


For the case of straight cracks, we use the Lobatto-Chebyshev method to 
solve these integral equations. Let 


r p = cos 


(? ~ I)* 


m. 


m — 1 


/or p = 1, 


(33) 
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The corresponding weights are 


7T 


W l = W m = 


2 (n - 1) 

The collocation points are 

(2 j ~ 1 )tt 


™r 


2(m — 1) 


for r = 1, m — 1 (34) 


fj — cos 


2m — 2 


/or / = 1, m — 1 


(35) 


For n cracks we get 2nm system of equations 
J 2 [ w pfv( T p_± + ki(T p ,{j,r Xl ,r Xj ,- • .)/ 2a (r p )u; p + k 2 (r p , tj,r Xl ,r X3 , • • ')f 22 (Tp)w p 

P = i 7r ( 7 > ~ 0) 

• • • + ' *)/ni(r p )wp + h(T,t,r Xl ,r Xn ,- • -^(t,,)™,,] = 4a n r^J 


/or j = 1, m — 1 


(36) 


^ \ 4" ^(^p) £j> r ai > > " ' ')f 2 l(' ! 'p) w p 4” ^-4 (^p) ) r *i > T’xj ) 

p=l 7r ( T P ~ 0) 


■)f 22 (r p )w p 


4~ ^3 ( T P) £77 r x I J ^Xn) ' ■ ')fnl{ T p} W p 4" ^4 (r p , £j’, Pa-J , P Xti , • • ”)/n2 ( T p)' u; p. — 40110”, 


IT 

W 


/or / = 1, m — 1 


(37) 


Z)[^i( r p>6» r *i> r * 3 .-- *)/2 i( t pM> 4” ^(7>,^,r Xl ,r Xl ,- • ■)fz 2 {T p )w p + 

P = i 7r l T p - 4i) 


4- ^i(7-,^,r xi ,r x „, •••)/„! (r p )m p -f k 2 (r,t,r Xl ,r Xn ,- - • )f n2 (r p )w p ] = 4 a n r : 


2T 

xy 


/or j = 1 , m — 1 


(38) 


r x i > r x 3 > * * ‘)f 2 i( T P ) w p 4~ Ar 4 (r p> r gl , r aa , • • •)f 22 (r‘p)w p -j- . 

x (. t p “ ?;) 


p=i 


4" ^(^pj^jj^xi j 2 "x n j ■ ' ‘)fnl(,' r P ) w p 4" ^(^pj £7 j r xi j r x n } " ‘ ')fn2{Tp)w p 40n0” ] 


2T 

yy 


/or / = 1, 


772 


(39) 



p-i 

+M T >£ 


X>[^ 3 ( t p’ 

P =1 

*f ^3 ( 7p ) 


0 > r x, ) T-x, i ' * *)/ll( T p) lu p "J* ^“2 (^pi ) r xi ) T’xj > * " ■)/l2( 7 p)' li, p 


> r x , , T Xn , ' 


•)/2i( t p)^p + k 2 (T p ,{j,r Xl ,r Xn ,- ■ ■)f 22 (r p )w p 


• + ] = 4a n r a 

7r ( r p-£>) 

/or j = 1, ..., m — 1 


>/nl (^p) 


nr 

*V 


'.ji T x n t r x u‘ ‘ ")f\l( T p) w p + k^Tp, £j, J* Xn , r Xj , ■ • ‘)f 22 (Tp)w p 

> ’’xni > ‘ ‘ *)/2l( 7 p)' li, p “t" ^4 (^p) /j> r Xn j 7-xj > * " ")f22(Tp)w p 

| ' U, p/n 2 ( r p)i _ 4 nT 

' ' “ 11 » 

for j = 1, .... m — 1 


/liOvK = 0 

p=i 


£ /l2(7>K = 0 

p=l 


E M t p) w p = 0 

P =1 


E M r rK = 0 

p=i 


E /mfoH = 0 

p= i 
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E M t p K = o 

p=i 


(47) 


We write the system of equations in matrix form 


[>!)(/) = (?) 


(48) 


Let 


JP — lj *i>j 777. 

J — lj >-«j 777. 1 

V P ~ *•*> m J 


_!£J_ 


*(n-6) 




\ wx 


mm 


\ 


7 r ( T m“^m-l ) 


) 


m*m 


and 


/ w p ki(r p ^j, \ 
r *«i r iu •••) 

P = l,-*™ 

V j = 1 rn - 1 


’’x.i r S(,> ••• ) ••• , I't,] •••) \ 


\ w lki{Th £m-l j j •••) ••• a, ra^i( 1 rn) ^m-li •■•) / (m — 1)» 


where i = 1, ...4; a, 6 6 1, ...77; and a ^ b. 


The matrix [A] is 



10 






















11 




12 



In order to calculate matrix [A], we calculate diagonal block matrix 




13 


and fill it into diagonal blocks in matrix [A]. 

Calculate the k{ block matrix and fill them into proper places in [A] for i = 1, 


ki(block ) = 


/ 


uip A:,' ( Tp , 

™pki{r p ,ti, 

™ p ki{T p ,£j, \ 


0 

r *l > ^**3 5 •") 

... r Xl , r Xn _ t , ...) 

r *l 3 r x„ ) •••) 



p = 1, m 

p = 1 , m 

p = 1 , m 



J = 1 

j = 1 , m — 1 

j = 1 , m — 1 


w pki(r p) (j, 

w p ki(r p , (j, 

w p ki{r p , £ji 

\ 


r *n) r il J •") 

r 2») r *l 1 "•) 

•" r *n> r *n-l 5 **•/ 


p = 1, m 

p = 

p = 1, ...,m 

0 

^ y = 1 

j = 1 

J — 1 , • , 771 “ * 1 

/ 


(nm— 2n)*nm 


After solving the system of equations, we get the stress intensity factors 


^’(l) = 1) («) 

4a u 

*F(i) = (so) 

4a u 

fcF’(-l) = — T~/j2(— 1) (51) 

4au 

#(-i) = -^-/»i(-i) (52) 

4a n 


n. 


where j = 1, 



FORTRAN PROGRAM FOR 


CALCULATING THE STRESSES OF n CRACKS 
EMBEDDED IN ISOTROPIC PLATE 

Main Program 

1) Call the subroutine indata to input all the informations of each crack 
and the far field stress. 

2) Select the collocation points. 

3) Call the subroutine lobatto to calculates the abscissae weights and 
collocation points by using Lobatto- Chebyshev integral technique. 

4) Call subroutine solsy to construct matrix A and right hand side of 
the system equations 

5) Call system subroution dlslrg to solve the large system of equations. 

6) Calculate the stress intensity factors. 

Subroutine indata 
An interface with the user. 

1) Tell user to input the elastic modulus E. 

2) Tell user to input the position ratio v. 

3) Let the user to choose the problem of plane stress or plane strain. 

4) Input the far field stress information. 

5) Do loop from 1 to n to input the position, angle and length of each crack. 



2 


6) Display all the parameters in the tables for user to make changes 
according to the index number. 

Plot all cracks proportionally on the screen. 

1) Initialize the crack representation. The ith crack is a line segment 
composed of number i. 

2) Calculate the maximum absolute value of the x— axis and y— axis involved 
in the real problem to decide the proportion of the graph which will be 
displayed according to the real problem. 

3) Do loop from 1 to n plot the center and tips of the ith crack 
and fill the middle with the number i for ith crack. 

Subroutine lobatto 

1) Calculation of the abscissae r. 

2) Calculation of the weights w. 

3) Calculation of the collocation points. 

Subroutine solsy 

1) Calculate diagonal blocks of block matrix A. 

2) Using the formulas of the four kernels derived from macsyma to calculate 
the (mn-2n)*nm block matrices, where n is the number of cracks, 

m is the collocation points. 

3) Fill block matrix A with the four block matrices in proper positions. 



3 


4) Calculate the right hand side of the system of equations by using the 
boundary conditions. 

5) Call the system subroutine dlslag to solve the system of equations. 

6) Calculate the stress intendity factors. 



* * ABT * * 

•k 

* Ten or less Cracks in the isotropic infinite plate 

k 


c 

c 


* 

* 


double precision gt (50) ,gc (50) , w (50) 

double precision f un (5000 ) , a (5000 , 5000 ) , p (50 00 ) 

double precision rx (10 ) , ry (10) , gf (10 ) , al (10) 

integer ncr 

double precision E,gn 

double precision gsOxx, gsOyy, gtOxy 

double precision all 

integer n 


common/property/ all 

common /numcrack/ ncr 

common /matprop/ E, gn 

common /load/ gsOxx, gsOyy, gtOxy 


external lslrg 
common /worksp/ rwksp 
real rwksp (9000000) 
call iwkin ( 9000000 ) 


The Problem of Less Than 11 Cracks 
in Infinite Isotropic Plate 


######' 

######' 


print *, ' #####################################################' 
print *,'###### 
print *,'###### 

print *, ' #####################################################' 
write (*, 104) 

************************************************************** 

* to change max number of cracks more than 10 

* the dimension of rx, ry, gf , al, . . . must be increased 
************************************************************** 

* 

* 


reading the problem data 


parameter (ncr=2) 

print Input the number of cracks less than 11' 
read (*, *) ncr 

print *, 'what, your input is ' , ncr 
write (*,104) 

call indata (rx, ry, gf, al) 


* V V w w w w w w vvw w w w w w w w w w vwvwwww w w w w w w w w w w 

* here we put an expert logic for the 

* optimum collocation points number 

^AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA 


write (*,*)' input number of collocation points n (less then 50)' 
read (*, *) n 

print *,'what your input is ',n 


* Here we put the choice of collocation technique (later) 

★ 

* calculating the collocation points, abbscissae and weigths 


* o 



* using Lobatto-Tchebyshev integral technique 
call lobatto (n, gt , gc, w) 

* 

write (*, 104) 

104 format (3x,/) 

* solve the system of equations and 

* calculate the stress intensity factor 

call solsy (n, rx, ry , gf , al, gt , gc, w) 

stop 

end 

k 


★ 

* 

k 

k 

k 

k 

k 

k 

k 

k 

k 


Subroutine for the problem setup 

- type of the problem (plain stress or plain strain) 

- number of cracks (ncrack=2) <this should be send from symb> 

- material properties (E, Nu=gn) 

- geometry of the crack locations (rx, ry , f i=gf ) 

- crack length (half of a crack=al) 

- far field loading condition in global coordinates 

0 0 0 

( Sigma = gsOxx Sigma = gsOyy Tau = gtOxy ) 
xx yy xy 


kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

subroutine indata (rx, ry, gf, al) 

character*l line (1 : 70, 1 : 30) 
character*l c(10),ans 

integer x,ncr,i, j, k, 1, kl, 11, 111, kkk, num 
double precision rx (ncr) , ry (ncr) , gf (ncr ) , al (ncr) , E, gn 
double precision gsOxx, gsOyy , gtOxy 
double precision all, xx, g (ncr) 

* 

double precision max,maxx,maxy,mxx,mxy 

k 

common/property/ all 
common/numcrack/ncr 
common/matprop/E, gn 
common/load/gsOxx, gsOyy, gtOxy 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

k 


k 

k 


Input data 


★ 


pi=dacos (-1 . OdO) 

k 

k 

k 

p j- j_nt * / r kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkf 

write (*,103) 

print *,' Input the elastic modulus E ' 
read(*, *) E 

print *,'what your input is ' ,E 
write (*, 103) 

print *, ' ******************************************************' 
write ( *, 103) 

k 




print *,' Input Possion ratio v ' 
read(*, *) gn 

print *,'what your input is ' , gn 
write ( *, 103) 

pj- 2_ nt * f r *********************************'******'**'*************' 

write { *, 103) 

print *, 'Choose the problem of plane stress or plane strain' 
print *,'For plane stress input 1' 

print *,'For plane strain input an integer other than 1' 
read(*, *) x 
IF (x .EQ. 1) THEN 
all=E 

print *,'what you select is plane stress problem' 

ELSE 

all=E/ (l-gn**2) 

print *,'what you select is plane strain problem' 

END IF 
write (*, 103) 

print * f r *********************************************'*'**' : * : '**'*■'fc'* , 
write (*,103) 

print *, ' ####################################################' 
print *,'##### Input the far field stress information #####' 
print *, ' ####################################################' 
write ( *, 103) 

pjfint * f r *'k********'k*'k*'kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkf 

write ( *, 103) 

print *,' Input the far field stress field gsOxx' 
read(*, *) gsOxx 

print *,'what your input is ', gsOxx 
write (*, 103) 

print * , ' ******************************************************' 
write (*, 103) 

print *,' Input the far field stress field gsOyy' 
read(*, *) gsOyy 

print *,'what your input is ', gsOyy 
write (*, 103) 

print * f r ************************************'*'' : * : ■*** : * r *' : * : '' : * ^ **■*'*'*■'*'*** , 
write ( *, 103) 

print *,' Input the far field stress field gtOxy' 
read(*, *) gtOxy 

print *,'what your input is ', gtOxy 
write (*, 103) 

k ^ t kk kkkkkkkkkk k k k k k k k k kk k k kk k k k k k k k k k k kk kk k k k k k k k k k k kkk ~k t 

write (*, 103) 


do 200 i=l,ncr 

print *,'###################################################' 
print *,'##### Input the information of the crack #',i,' ##' 

print *, ' ###################################################' 
write { *, 103) 

pj-^j-j-^- -k ^ t kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk' 

write ( *, 103) 

print *,' Input the crack position vector rx of crack #',i 
read ( *, *) rx (i) 

print *,'what your input is ',rx(i) 
write (*, 103) 

p2-j_j-j^- k f r kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk' 

write (*, 103) 

* 



print Input the crack position vector ry of crack #',i 
rea d(*, *) ry (i) 

print *,'what your input is ' ,ry(i) 
write ( *, 103) 

pj- 2 _ nt * , 9 ******************************************************* 

write (*, 103) 

print Input the crack rotation gf of crack #',i 
read(*, *) gf (i) 

print *,'what your input is ' ,gf(i) 
g(i)=gf (i) 

gf (i)=(2*gf (i)*pi)/360 
gf ( 1) =dacos (dcos (gf (1) ) ) 
print *,'In arc gf (i) =' , gf (i) 
write (*, 103) 

p f r ******************************************************* 

write ( *, 103) 

print Input half of the length of the crack #',i 
read(*, *) al (i) 

print *,'what your input is ',31(1) 
write ( *, 103) 

★ * * ******************************************************* 
write ( *, 103) 

200 continue 

103 format (3x,/) 


★ ^•★★*****]_o********20*** : *'**'**30 ********40 ********50 ******** 60 ********70* 

★ 

* Print out a table of input data 

* 

*********10********20********30********40********50********60********70* 

* 


400 


print 
print 
print 
print 
print 
print 
print 
print 
print 
print 
IF (x 
print 
ELSE 
print 
END IF 

print * , ' * * ' 

"k f f kkkk'kickkkkkiirkkkkkkkkkkkkk'kkkkkkkkkkkkkkkkkkkkklekkkkkkkr 

write ( *, 103) 

print *,'Do you want to change any of the value above? ' 
print Answer y for yes, n for no. ' 
read (*, *) ans 

IF (ans .EQ. f y' ) THEN 

print *,' Input the number to represent the element' 
print *,'you want to change. ' 
read (*, *) num 
IF (num .EQ. 1 ) THEN 

print Input the elastic modulus E ' 
read (*, *) E 

print *,'what your input is ',E 
write (*, 103) 

END IF 

IF (num .EQ. 2) THEN 


* t * * *' 

*, ' * 1 * E * ' , E 

* r * * *r 

r 

■k f **************** ********* ****** **** ************* 

* f ' * * *' 

*, 9 * 2 * v * ' , gn 


★ r * 


* r 


f f ******************************************************' 

* r * * r 

r 

.EQ. 1) THEN 

*,'* 3 * Problem of plane stress' 


*,'* 3 * Problem of plane strain' 



401 


print ' Input Possion ratio v ' 
read(*, *) gn 

print *,'what your input is ' , gn 
write (*, 103) 

END IF 

IF (num .EQ. 3) THEN 

print *, 'Choose the problem of plane stress or plane strain' 
print *,'For plane stress input 1' 

print *,'For plane strain input an integer other than 1' 
read {*, *) x 
IF (x .EQ. 1) THEN 
all=E 

print *,'what you select is plane stress problem' 

ELSE 

all-E/ (l-gn**2) 

print *,'what you select is plane strain problem' 

END IF 
write (*, 103) 

END IF 
goto 400 
END IF 
write (*, 103) 

pj- jQ-J- k ^ t kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkr 


print 

k f k ★ 

t 

0 

* r 


print 

k t -k 4 k 

Sigma 

k / 

r 

gsOxx 

print 

k t * k 

t 

XX 

kf 


print 

* f 9 ******************************************************' 

print 

k r k k 

/ 

0 

kr 


print 

k f k 5 * 

Sigma 

k r 

r 

gsOyy 

print 

k f k k 

r 

yy 

k r 


print 

k ^ f kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkf 

print 

k r k k 

f 

0 

kr 


print 

*,' * 6 * 

Tau 

k r 

r 

gtOxy 

print 

k r k k 

t 

yy 

k r 


print 

k f f kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk t 

write (* 

, 103 ) 




print 

*,'Do you want to change any of the value above? ' 

print 

*, ' Answer y 

for yes, n 

for 

no . ' 

read ( * 

, *) ans 




IF 

(ans .EQ. 'y 

' ) THEN 




print *,' Input the number to represent the element' 
print *,'you want to change. ' 
read (*, *) num 
IF (num .EQ. 4) THEN 

print Input the far field stress field gsOxx' 
read(*, *) gsOxx 

print *,'what your input is gsOxx 
write (*, 103) 

END IF 

IF (num .EQ. 5) THEN 

print Input the far field stress field gsOyy' 
read (*, *) gsOyy 

print *,'what your input is gsOyy 
write (*, 103) 

END IF 

IF (num .EQ. 6) THEN 

print Input the far field stress field gtOxy' 
read (*, *) gtOxy 

print *,'what your input is ', gtOxy 
write (*, 103) 

END IF 
goto 401 
END IF 

pause '******* Hit the key <return> to continue 
write (*, 103) 
do 201 i=l,ncr 


kkkkkkk r 



402 


j=ll+4* (i-1) 

print * f r ******************************************************' 
print * , ' * * * ' 

print *,'* ',j,' * rx ( ' , i, ' ) * rx(i) 

print *, ' * * *' 

print * f 9 ■k-k-k-k'k-k'k-k-k-k’k-k-k'k-k-k-k-k-k-k-k-k-k-k'k-k'kic-k'k-k-k'k'k'k-k'k'k'k'k'k'k'k'k'k'k-k-k-k'k'k-kic-kf 

print * , ' * * *' 

print *,'* ',j+l,' * ry(',i,' ) * ', ry(i) 

print *, ' * * *' 

print *, ' * * *' 

print *,'* ' , j+2 , ' * gf(',i,' ) * ', g(i) 

print * , ' * * *' 

print 

print *,'* * *' 

print *,'* ' , j+3, ' * al(',i,' ) * ', al(i) 

print *, ' * * *' 

f t 'k-k'k-k'k-k'k-k-k-kicic'k'k'k'k-k'k-k-k'k'k'k'k'k'k'k'k'k'k'k'k'k-k-k'k-k-kick-k'k-k'k'k-k-kic'k'k-k'k'k'kf 

write (*, 103) 

print *,'Do you want to change any of the value above? ' 
print Answer y for yes, n for no. ' 
read(*, *) ans 

IF (ans .EQ. ' y' ) THEN 

print *, 'Input the number to represent the element' 
print *,'you want to change. ' 
read(*, *) num 
IF (num .EQ. j ) THEN 

print Input the #',i,' crack position vector rx(',i,' )' 

read (*, *) rx (i) 

print *,'what your input is ',rx(i) 
write (*, 103) 

END IF 

IF (num .EQ. j+1 ) THEN 

print Input the #',i,' crack position vector ry(',i,' )' 

read (*, *) ry (i) 

print *,'what your input is ',ry(i) 
write (*, 103) 

END IF 

IF (num .EQ. j+2) THEN 

print *,' Input the #',i,' crack rotation gf(',i,' )' 

read (*, *) gf (i) 

print *,'what your input is ',gf(i) 
g(i)=gf (i) 

gf (i) = (2*gf (i)*pi)/360 
print *,'In arc gf ( i) =' , gf (i) 
write (*, 103) 

END IF 

IF (num .EQ. j+3) THEN 

print Input half of the length of the #',i,' crack' 
read (*, *) al (i) 

print *,'what your input is ',al(i) 
write ( *, 103) 

END IF 
goto 402 
END IF 
201 continue 

write (*, 103) 

write (*, 103) 
write (*, 103) 

* 

* 

* 

*********10********20********30*******M0********50********60********70* 

★ 

* Initialize the screen with all blank 



■k 

do 3 i = 1,70 
do 3 j=l,30 
lined, j) = ' ' 

3 continue 

k 

* Create a coordinate on the screen 

* 

do 5 i=l, 70 

line (i, 15) = ' 

5 continue 
line (5,15)=' | ' 
line (15,15 )=' [' 
line (25, 15)=' | ' 
line (35, 15)=' | ' 
line (45, 15)=' | ' 
line (55,15)=' | ' 
line (65, 15)=' | ' 

do 6 i=l,30 

line (35, i) =' I' 

6 continue 

line (35,5)='-' 
line (35,10)=' 
line (35, 15) =' | ' 
line (35,20)='-' 
line (35,25)='-' 

* 

*********20********20********30********40********50********60********70* 

★ 

* Initialize the crack representation 

★ 

c(l)='l' 
c (2 ) =' 2 ' 
c (3 ) =' 3 ' 
c (4 ) =' 4 ' 
c (5) =' 5' 
c (6)=' 6' 
c (7 ) =' 7 ' 
c (8 ) =' 8 ' 
c (9)=' 9' 
c (10) =' $' 

★ 

k 

* Calculate the maximum value of x-axis 

* and y-axis involved in the real problem 

k 

*********io********20********30*******MQ********50********60********70* 

* 

max=0 

do 8 i=l,ncr 

if ( abs (rx (i) ) +al (i) .GT. max ) then 
max=abs (rx (i) ) +al (i) 
end if 

if ( abs (ry (i) ) +al (i) .GT. max ) then 
max=abs (ry (i) ) +al (i) 
end if 

8 continue 

k 

kkkkkkkkk IQ kkkkkkkk 20 ******** 3()kkkkkkkk 4Q k*kkkkkk $0 kkkkkkkkQQkkkkkkkklQk 
k 
k 


Calculate the crack center on the screen 





*********20********20********30********40********50********60********70* 
do 10 i=l,ncr 
k=INT (rx (i) * (30/max) ) 

* print *, ' k=' , k 
1=INT (30*ry (i) /max) / 2 
line (35+k, 15-1) =c(i) 

*********io********20********30********40********50********60********70* 

★ 

* Calculate right end of the crack on the screen 

* 

*********3_Q******-*r*20********30********40********50********60********70* 
k=INT ( (rx (i) +dcos (gf (i) ) *al (i) ) *30/max) 

1=INT (30* (ry (i) +dsin (gf (i) ) *al(i) ) /max/2) 
line (35+k, 15-1) =c(i) 

* 

*********]_o********20********30********40********50********60********70* 

~k 

* Plot the line between crack center and the right end (x-axis) 

* 

*********]_q********20********30********40********50 ******** 60********70* 

★ 

if (rx (i) *30/max .LT. k ) then 
kkk=l 

kl=rx (i) * (30/max) +1 
do while (kl .LT. k) 

11= (30*ry (i) /max+kkk*dtan (gf ( i ) ) ) / 2 
line (35+kl, 15-11) =c (i) 
kl-kl+1 
kkk=kkk+l 
end do 
else 

kkk=l 

kl=rx (i) * (30/max) -1 
do while (kl .GT. k) 

* ll=-14*ry (i) /maxy- (kl-rx (i) * (34/maxy) ) *14/34*dtan (gf (i) ) 

* 11=14* (ry (i) -kkk*maxx/34*dtan (gf (i) ) ) /2/maxy 
11= (30*ry (i) /max-kkk*dtan (gf ( i ) ) ) / 2 

line (35+kl, 15-11) =c (i) 
kl-kl-1 
kkk=kkk+l 
end do 
end if 

* 

*********j_0********20********30********40********50********60********70* 

★ 

* Plot the line between crack center and the right end (y-axis) 

★ 

*********]_0********2O********3O********4O********5O********6O********7O* 

★ 

if (ry (i) *30/max/2 .LT. 1 ) then 
kkk=2 

ll=ry (i) *30/max/2+l 
do while (11 .LT. 1) 

kl=30*rx (i) /max+kkk/dtan (gf (i) ) 
line (35+kl, 15-11) =c(i) 

11 = 11+1 
kkk=kkk+2 
end do 
else 
kkk=2 

ll=ry (i) *30/max/2-l 
do while (11 .GT. 1) 

kl=30*rx (i) /max-kkk/dtan (gf (i) ) 
line (35+kl, 15-11) =c(i) 

11 = 11-1 



kkk=kkk+2 
end do 
end if 


*********iq********20********30********40********50********60********70* 

k 

* Calculate the left end of the crack on the screen 

★ 

*********^0********2Q***'*'*'*'**3O********4O********5O********6O********7O* 

★ 

k=INT ( (rx (i) -dcos (gf (i) ) *al (i) ) *30/max) 

1=INT (30* (ry (i) -dsin (gf (i) )*al(i) )/max/2) 
line (35+k,15-l)=c(i) 

★ 

*********]_q********20********30********40********50********60********70* 

k 

* Plot the line between crack center and the left end (y-axis) 

★ 

*********10********20********30********40********50********60********70* 

★ 

if ( 1 .LT. ry (i) *30/max/2 ) then 
kkk=2 

ll=ry (i) *30/max/2-l 
do while (1 .LT. 11) 

kl=30*rx (i) /max-kkk/tan (gf (i) ) 
line (35+kl, 15-11) =c (i) 
kkk=kkk+2 
11 = 11-1 
end do 
else 
kkk=2 

ll=ry (i) *30/max/2+l 
do while (1 .GT. 11) 

kl=30*rx (i) /max+kkk/tan (gf (i) ) 
line (35+kl, 15-11 )=c(i) 
kkk=kkk+2 
11 = 11+1 
end do 
end if 

* 

*********]]_q;A:*******20********30********40********50********60********70* 

★ 

* Plot the line between crack center and the left end (x-axis) 

* 

*********]_o********20********30********40********50********60********70* 

★ 

if ( rx(i)*30/max .LT. k ) then 
kkk=l 

kl=rx (i) * (30/max) +1 
do while (kl .LT. k) 

11= (ry (i) *30/max+kkk*dtan (gf (i) ) ) /2 
line (35+kl, 15-11 )=c (i) 
kl=kl+l 
kkk=kkk+l 
end do 
else 
kkk=l 

kl=rx (i) * (30/max) -1 
do while (kl .GT. k) 

11= (ry (i) *30/max-kkk*dtan (gf (i) ) ) /2 
line (35+kl, 15-11 )=c(i) 
kl=kl-l 
kkk=kkk+l 
end do 
end if 

* 



-k 

10 continue 

k 

k 

*********]^********20********30********40********50********60********70* 

k 

* Show the graph on the screen 

★ 

*********io********20********30*******MO********50********60********70* 

k 

k 

do 4 j=l, 30 

print *, (line (i, j) , i = 1,70) 

4 continue 

* 

*********io********20********30********40********50********60********70* 

* print *, 'rx(l) = ' ,rx (1) , ' ry (1) = ' , ry (1) 

* print *, 'rx (2) = ' , rx (2) , ' ry (2) =' , ry (2) 

* print *, 'gf (1) =' , tan (gf (1) ) , ' gf (2) =' , tan (gf (2) ) 

* print *, ' al (1) =' , al (1) , ' al (2) = ' , al (2) 

*********i0********20********30*******M0********50********60********70* 

* 

* the end of the subroutine of input data 

★ 

* 11 format (lx, 61al, ')' ) 

write (*, 103) 

pause '******* Hit the key <return> to continue *******' 
write (*, 103) 
return 
end 

* 



* 


* Subroutine for collocation points, abscissaes, and weights 

* using Lobatto-Tchebyshev integral formula 

* 



subroutine lobatto (n, gt , gc, w) 
integer n 

* 


double precision gt (n) , gc (n-1) , w (n) , pi 

* 


* 

* Calculation of the tau - abscissae 

★ 

pi=dacos (-1 . OdO) 

* print *, ' pi' , pi 

* write (*,105) 

do 10 j=l,n 

gt ( j) =dcos ( ( ( j — 1 ) *pi) / (n-1) ) 

10 continue 

* 

* Calculation of the weights 

k 

w (1) =pi/ (2* (n-1) ) 
w (n) =w (1) 
do 20 i=2,n-l 





w (i) =pi/ (n-1) 

20 continue 

* 

* Calculation point x 

* 

do 30 j=l , n-1 

gc ( j)=dcos ( ( (2* j— 1 ) *pi) / (2*n-2) ) 

* print *,'gc',gc(j) 

30 continue 

★ 


105 format (3x,/) 

★ 


return 

end 

************************************************************************ 

★ 

* Subroutine for solving the system of integral equations 

* 

*jc*Jr*ic1ric*jc1cJc*icir**jr*icicJc*icie*jrjr1cJciejcjrJcicicicJc*1eirie'k'k'k'k-k-k'k-kic-k'k'kic-k-kicicicjcjcjrjc1cicic1c1c1c1tir 

★ ★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★★■A-***** 

subroutine solsy (n,rx, ry, gf , al, gt , gc, w) 

★ 

integer i, j , k, 1, n, ncr 

double precision krl (ncr*n, ncr*n) , kr2 (ncr*n, ncr*n) 

double precision kr3 (ncr*n, ncr*n) , kr4 (ncr *n, ncr *n) 

double precision rx (ncr) , ry (ncr) , gf (ncr) , al (ncr) , gt (n) , gc (n-1) 

double precision w (n) , E, gn, gsOxx, gsOyy, gtOxy 

double precision fun (2*ncr*n) ,a (2*ncr*n, 2*ncr*n) ,p (2*ncr*n) 

double precision kilpl (ncr) , kilml (ncr ) , ki2pl (ncr) , ki2ml (ncr ) 

double precision pi, ww, wwm, parti, part2 

double precision gh,plala2, p2a2 

double precision pl,p2 

double precision gsyiyi (ncr ) , gtxiyi (ncr ) 
double precision pixi (ncr) , qixi (ncr) 
double precision all 

* 

double precision rlx, rly, r2x, r2y, gf 1, gf 2, al, a2 

★ 


★ 

★ 

* 

* 


common/property/ all 
common /numcrack/ ncr 
common/matprop/E, gn 
common/load/gsOxx, gsOyy, gtOxy 

pi=dacos (-1 . OdO) 


print gsOxx' , gsOxx 

print gsOyy' , gsOyy 

print ' gtOxy' , gtOxy 
print *, 'E' ,E 
print * , ' v' , gn 
print *, 'ncr' ,ncr 
print *, ' rx (1) ' , rx (1) 
print ' rx (2) ' ,rx (2) 
print *, ' n' , n 
print *,'gt',gt(l) 



* 

do 66 

j=l,n-l 

* 

print 

' gc' , gc ( j) 

* 

print 

*, 'all' ,all 

* 66 

continue 

A: 

★ 

print 

*,'w',w(l) 



★ 


do 292 i=l,ncr*n 
do 293 j=l,ncr*n 
krl (i, j) =0 
293 continue 

292 continue 

* 


* 

* Initialize the matrix [kr2] 

* 


* 

do 294 i=l,ncr*n 
do 295 j=l,ncr*n 
kr2 (i, j ) =0 
295 continue 

294 continue 

* 

************************************************************************ 
' * 

* Initialize the matrix [kr3] 

* 


* 

do 296 i=l,ncr*n 
do 297 j=l,ncr*n 
kr3 (i, j) =0 
297 continue 

296 continue 

* 


* 

* Initialize the matrix [kr4] 

* 


★ 

do 298 i=l,ncr*n 
do 299 j=l,ncr*n 
kr4 (i, j) =0 
299 continue 

298 continue 

* 


★ 

* Initialize the matrix [A] and [P] 

* 


* 


do 300 i=l,2*ncr*n 







301 


do 301 j=l,2*ncr*n 
a (i, j)=0 
continue 
P (i) =0 
300 continue 


★ 

★ 

* Calculate the diagonal elements of matrix [A] 

★ 

★ 


do 305 k=0, (2*ncr-l) 
do 40 i=l, n-1 
do 50 j=l,n 

a (k*n+i, k*n+j) =w( j) / (pi* (gt ( j) -gc (i) ) ) 

* 

* print *, ' adiagonal=' , a (k*n+i, k*n+ j ) 


50 continue 
40 continue 
305 continue 


* 

k 


do 306 k=0, (2*ncr-l) 
do 60 i=l,n 

a ( (k+1) *n, k*n+i) =w (i) 

* print ' adiagdown=' , a ( (k+1) *n, k*n+i) 

60 continue 

306 continue 

* 

★ 

★ 

★ 

********* krl ***** (3,1) ***** k2111 ********** 

* 

* 


do 310 k=l,ncr 
do 311 1=1, ncr 

* 

if (k .NE. 1) then 

* 

r2x=rx (k) 
r2y=ry (k) 
gf 2=gf (k) 
a2=al (k) 

* * 

rlx=rx (1) 
rly=ry (1) 
gfl=gf (1) 
al=al (1) 

* 

gh=gf 2-gf 1 

* 

pl= (r2y-rly) *dsin (gf 1) + (r2x-rlx) *dcos (gf 1) 

* 

p2= (r2y-rly) *dcos (gfl) - (r2x-rlx) *dsin (gf 1) 

* 

do 78 i=l , n-1 
do 79 j=l,n 

* 

plala2 = -pl+al*gt ( j) -a2*dcos (gh) *gc (i) 
p2a2 = p2+a2*dsin (gh) *gc (i) 
ww = p2a2**2+plala2**2 



wwm = plala2**2-p2a2**2 

parti = 2*dcos (2*gh) *plala2*p2a2**2-dsin (2*gh) *p2a2*wwm 
parti = partl/ww**2 

part2 = (-dsin (2*gh) *p2a2-dcos (2*gh) *plala2) /ww 

★ 

krl { (k-1) *n+i, (l-l)*n+j)=(al/pi)* (partl+part2 ) *w ( j ) 

* print *, ' krl=' , krl ( (k-1) *n+i, (1-1) *n+ j) 

* 

79 continue 

78 continue 

* 

end if 

* 

311 continue 
310 continue 

k 

k 

kkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkkk 

k 

kkkkkkkk kr2 ****** (3,2) ***** k2112 ********** 

* 

************************************************************************* 

* 

do 312 k=l,ncr 
do 313 1=1, ncr 

* 

if (k .NE. 1) then 

* 

r2x=rx (k) 
r2y=ry (k) 
gf2=gf (k) 
a2=al(k) 

** 

rlx=rx (1) 
rly=ry (1) 
gfl=gf (1) 
al=al (1) 

* 

gh=gf 2-gf 1 

* 

pl= (r2y-rly) *dsin (gfl) + (r2x-rlx) *dcos (gfl) 

* 

p2= (r2y-rly) *dcos (gfl) - (r2x-rlx) *dsin (gfl) 

* 

do 80 i=l,n-l 
do 81 j=l,n 

* 

plala2 = -pl+al*gt ( j ) -a2*dcos (gh) *gc (i) 
p2a2 = p2+a2*dsin (gh) *gc (i) 
ww = p2a2**2+plala2**2 
wwm = plala2**2-p2a2**2 

parti = dcos (2*gh) *p2a2*wwm+2*dsin (2*gh) *plala2*p2a2**2 
parti = partl/ww**2 

* 

kr2 ( (k-1) *n+i, (1-1) *n+ j ) = (al/pi) *part l*w ( j) 

* print *, ' kr2 = ' , kr2 ( (k-1) *n+i, (1-1) *n+ j ) 

* 

81 continue 

80 continue 

* 

end if 

* 

313 continue 

312 continue 


* 










if (k .NE. 1) then 


* 

do 356 i=l,n-l 
do 357 j=l,n 

★ 

a (2*k*n+i, 2*l*n+n+ j ) =kr2 (k*n+i,l*n+j) 

★ 

357 continue 
356 continue 

* 

end if 

* 

355 continue 
354 continue 

* 

* 

******** Input kr3 to [A] ************************************* 

* 

do 358 k=0,ncr-l 
do 359 l=0,ncr-l 

* 

if (k .NE. 1) then 

* 

do 360 i=l,n-l 
do 361 j=l,n 

* 

a (2*k*n+n+i, 2*l*n+j) =kr3 (k*n+i, l*n+ j) 

■Ar 

361 continue 
360 continue 

* 

end if 

* 

359 continue 

358 continue 

★ 

* 

******** Input kr 4 to [A] ★★★★★★*★■*•**★★★★★*****★■*★★★★★★**★*★★■*■■*• 

★ 

do 362 k=0,ncr-l 
do 363 l=0,ncr-l 

* 

if (k .NE. 1) then 

* 

do 364 i=l,n-l 
do 365 j=l,n 

* 

a (2*k*n+n+i, 2*l*n+n+ j ) =kr 4 (k*n+i, l*n+ j ) 

* 

365 continue 
364 continue 

★ 

end if 

★ 

363 continue 

362 continue 

* 


* 



* 


* 


* 


print *, 'matrix a= 
do 200 j=l,4*n,6 







print 502, ki2ml (k+1) , ki2pl (k+1) 

502 format (lx, ' k2 (-1) =' , el2 . 6, 3x, ' k2 (1) = ' , el2 . 6) 

★ 

391 continue 

★ 

★ 

return 

end 



/* solve the ODE get r */ 
f f : cc*%e A (r*y*s) $ 

fq: factor (dif f (ff*%e A (-%i*s*x),x,4) +2*dif f (ff*%e A (-%i*s*x),x,2,y,2) +dif f (ff*%e A (-%i- 

/* solve the characteristic equation */ 
ql : solve (fq*%e A (-r*y*s+%i*s*x) , ' r) $ 

/* it should be two identical real roots */ 
if length (ql) #2 then 'errorl$ 
if rhs (ql [1] ) #-rhs (ql [2] ) then 'error2$ 
r : abs (rhs (ql [ 1] ) ) $ 

/* the solution of the ODE at y>0 */ 

fxyp: (ccl+cc2*y) *%e A (-ev(r) *abs (s)*y-%i*s*x)$ 

/* the solution of the ODE at y<0 */ 

fxym: (cc3+cc4*y)*%e A (ev(r)*abs(s)*y-%i*s*x)$ 

/* calculate stresses gs[i] */ 

block (gs [1] : dif f (fxyp, y, 2 ) , gs [2] : dif f (fxyp, x, 2) , gs [3] : diff (fxyp, x, 1, y, 1) , 
gs [ 4 ] : diff (fxym, y, 2 ) , gs [5 ] : dif f (fxym, x, 2) , gs [ 6] : diff (fxym, x, 1, y , 1 ) ) $ 

/* the shear tractions q(x)=0 so gtxy=0 for y=0 . */ 

/* According to continuity conditions cc3 and cc4 could be expressed by ccl */ 
pp:diff ( (ccl+cc2*y) *%e A (-ev (r)*abs(s)*y)-(cc3+cc4*y)*%e A (ev(r)*abs(s)*y),y,l)$ 
ppp : rat subs t (0, y, ( (ccl+cc2*y) *%e A (-r*abs(s)*y)-(cc3+cc4*y)*%e A (r*abs(s)*y)))$ 

11s : solve (ppp, cc3) $ 

cc3 : rhs (11s [ 1 ] ) $ 

ppt : rat subs t (0, ' y,pp) $ 

11s : solve (ppt, ' cc4) $ 
cc4 :ev(rhs (11s [1] ) ) $ 

/* calculate auxilliary function fl (x) */ 
eqql : dif f ( (integrate (fxyp, x) ) , y, 2 ) $ 
eqq2 : diff (fxyp, x, 1) $ 

flxp-.diff ( (ratsubst (0,y, (all*eqql+al2*eqq2) ) ) ,x,l)$ 
eqml : dif f ( (integrate ( ' ' fxym, x) ) , y, 2 ) $ 
eqm2 : dif f ( ' ' fxym, x, 1) $ 

f lxm: dif f ( (ratsubst (0 , y, (all*eqml+al2*eqm2) ) ) , x, 1) $ 
mi : (ratsimp ( (f lxp-f lxm) / (2*%pi) } ) *2*%pi*%e A (%i*s*x) $ 
lx: solve (mi-fl (t) *%e A (%i*s*t) , cc2) $ 
cc2 : rhs (lx [ 1] ) $ 

/* calculate auxilliary function f 2 (x) */ 
eqql : diff ( (integrate (fxyp, y) ) , x, 2} $ 
eqq2 : diff (fxyp, y, 1) $ 

f2xp: ratsubst (s*abs (s) ,ks, (ratsubst (ks*abs(s),s A 3, (diff ( (ratsubst (0,y, (all*eqql+al2' 
eqml : dif f ( (integrate ( ' ' fxym, y) ) , x, 2 ) $ 
eqm2 : diff ( ' ' fxym, y, 1) $ 

fm2 : ratsubst (0, y, (all*eqml+al2*eqm2 ) ) $ 

f2xm: ratsubst (s*abs (s) , ks, (ratsubst (ks*abs(s),s A 3, (diff (fm2, x, 1) ) ) ) ) $ 
f2x : ' integrate ( (ratsimp ( (f2xp-f2xm) / (2*%pi) ) ) , s,minf , inf) $ 

/* solve ccl expressed by the auxiliary function */ 

mi : (ratsubst (s*abs (s) , ' integrate (s*abs(s)*%e A (-%i*s*x),s,minf,inf),f2x))*2*%pi$ 
lx : solve (mi=' integrate ( f 2 (t) *%e A (%i*s*t) ,t,-a,a) ,ccl) $ 
ccl : rhs (lx [ 1] ) $ 

ccl : ratsubst (f2 (t) *%e A (%i*s*t) , ' integrate (%e A (%i*s*t)*f2(t) ,t,-a,a) , ccl) $ 

/* calculate cc2 cc3 and cc4 expressed by the auxiliary function */ 
cc2 : expand (cc2 : xthru (ev (cc2 ) ) ) $ 
cc21 : last (cc2 ) $ 
cc2 : first (cc2 ) $ 

cc21 : ratsubst (abs (s ) , kb, (ratsubst (1, abs (s ) , (ratsubst (kb,s A 2,cc21) ) ) ) ) $ 
block (cc2 : cc2+cc21, cc3 : ev (cc3) , cc4 : ev (cc4 ) ) $ 

/* Get the answer of gsxxp expressed by the auxiliary function */ 
for j : 1 thru 3 do (mwp : expand (ev (gs ( j ])) , 



/* calculate |s| when s>0 */ 
mip : expand ( subs t ( s , abs ( s ) , mwp ) ) , 

/* calculate |s| when s<0 */ 

mim: expand (rat subst (-s, s, (subst (-s, abs ( s) ,mwp) ) ) ) , 

/* combine the two terms */ 

mi : ratsubst (ywm, s*wm, (ratsubst (ywp, s*wp, (subst (wm, %e A (-s*y+%i*s*x-%i*s*t) , (subst (wp 
if j=3 then ( mi : ratsimp (mi) ) , 

mi : distrib (subst ( z, (t-x) , ( (subst ( (-1/ (-y-%i* (t-x) ) ) , wm, (subst ( (-1/ (-y+%i* (t-x) ) ) , wp 

ki : factor (expand (xt hr u ( first (mi) +f irst (rest (mi, 2 ) ) ) ) ) , 
k j : factor (expand (xt hr u (first (rest (mi, 1) ) +last (mi) ) ) ) , 
mi: subst ( (t-x) , z, ki+k j) , 

/* Get the answer of stresses gs[i] expressed by the auxiliary function */ 
gs t j] : ' integrate (mi, t, -a, a) , 

LD I SP LAY ( gs [ j ] ) , 


ki : ki*4 *%pi*all , 
k j : k j *4 *%pi*all , 
mi : ki+k j, 

gsxxpf [ j] ; subst ( (t-x) , z,mi) 

) $ 

/★**** calculate four kernels of two cracks ***********************************/ 

j 

for i:l thru 3 do ( 

gsxxpf It 1 [i] : subst (tl, t , (subst (xl, x, (subst (yl, y, (subst ( 1, f 1 (t ) , (subst (0 , f 2 (t ) , (gsxx] 
gsxxpf lt2 [i] : subst (t2 , t , (subst (x2 , x, ( subst (y2 , y, (subst ( 1, f 1 (t ) , (subst (0 , f 2 (t ) , gsxxp: 
gsxxpf 2tl [i] : subst (tl, t, (subst (xl, x, (subst (yl, y, (subst (1, f2 (t) , (subst (0, f 1 (t) , gsxxp: 
gsxxpf 2t 2 [i] : subst (t2 , t , (subst (x2, x, (subst (y2 , y , (subst (1, f 2 (t ) , (subst (0 , f 1 (t ) , gsxxp: 
) $ 

kill (wwm) $ 

gsx2 : subst (a2*gt, t2, gsxxpf lt2 [1] ) $ 
gsy2 : subst (a2*gt, t2, gsxxpflt2 [2] ) $ 
gtxy2 : subst (a2*gt, t2, gsxxpf It 2 [3] ) $ 

gm21b [1] : - ( ' ' gsx2-' ' gsy2) /2*sin (2* (-gh) ) +' ' gtxy2*cos (2* (-gh) ) $ 

gm21b [3 ] : ( ' ' gsx2+' ' gsy2 ) /2- ( ' ' gsx2-' ' gsy2 ) /2*cos (2* (-gh ) ) -' ' gtxy2 *sin (2 * ( -gh) ) $ 

gsx2 : subst (a2*gt , t2 , gsxxpf 2t2 [1] ) $ 
gsy2 : subst (a2*gt, t2, gsxxpf 2 t2 [2] ) $ 
gtxy2 : subst (a2*gt, t2, gsxxpf 2t 2 [ 3] ) $ 

gm21b [2] :-( ,, gsx2- ,, gsy2)/2*sin(2* (-gh ) ) +' ' gtxy2*cos (2* (-gh) ) $ 

gm21b [4 ] : ( ' ' gsx2+' ' gsy2 ) /2- ( ' ' gsx2-' ' gsy2 ) /2*cos (2* (-gh) ) -' ' gtxy2*sin (2* (-gh) ) $ 

cs.lcul 3 .t 0 the four kernels 
for i:l thru 4 do ( 

gm21[i] : subst (wwm, p3a2al , '2-p4al A 2 , (subst (ww, p4al A 2+p3a2al A 2 , (subst (p4al , p4-al*gc*sii 
if i=l then ( 

gm21 [ i] : combine (mult thru (first (gm21 [i] ) ) +multthru (rest (gm21 [i] , 1) ) ) 

\ ' I 

xf a=2 then ( gm21 [i] : combine (gm21 [i] )) , 
if i=3 then ( 

n :combine ( (last (gm21 [i] ) ) + (multthru (first (gm21 [i] ) ) ) + (multthru (first (rest (gm21 [i] , 1 

ttl : first (n) , 

ts2 :part (last (n) , 1) , 

1 1 2 : { (trigs imp (trigexpand (factor (first(ts2)+first(rest(ts2, 1) ) ) ) ) ) +last (ts2 ) ) /ww, 
gm21 [i] : ttl+tt2 
) , 

if i=4 then (gm21 [i] : combine (gm21 [i] )) , 


/* display the four kernels k[i] after simplification ********************/ 
k [i] : (a2/ (4*%pi*all) ) *gm21 [i] , 



LDI SPLAY (k [ i] ) , 

/* where 

p3a2al=~p3+a2*gt-al*gc*cos (gh) 
p4al=p4~al*gc*sin (gh) 
ww=p4 a 1 A 2 +p3a 2 a 1 A 2 
wwm=p3a2al A 2-p4al A 2 
*/ 

kill (ww, wwm) , 

/*** generate FORTRAN code for calculating numerically **********************/ 

fortran (p3a2al : -p3+a2*gt { j ) -al*gc (i) *dcos (gh) ) , 

fortran (p4al :p4-al*gc (i) *dsin (gh) ) , 

fortran (ww:' p4al A 2+' p3a2al A 2) , 

fortran (wwm: ' p3a2al A 2-'p4al A 2) , 

fortran (parti : (subst (dsin, sin, (subst (dcos, cos, first (first (gm21 [i] ))))))) , 
fortran (parti : 'parti/ ( f ww A 2)), 

fortran (part2 : (subst (dsin, sin, (subst (dcos, cos, last (gm21 [i] )))))) , 
kill(nl,n2,n3, n, p3a2al, p4al, ww, wwm) 

)$ 



(d3) 


totallp 


(c4) /* solve the ODE get r */ 
ff: 'cc*%e A ( ' r*' y* ' s ) $ 

(c5) fq:factor( 

diff (ff *%e A (— %i*' s*' x) , 'x, 4) + 

2*dif f (f f *%e A (-%i*' s* ' x) , ' x, 2 , ' y, 2 ) + 
dif f (ff *%e A (-%!*' s*'x) , 'y, 4) 

) ; 

2 24 rsy-%isx 

(d5) cc (r - 1) (r + 1) s %e 

(c6) remvalue (f f ) $ 

(c7) /* solve the characteristic equation */ 
ql : solve (fq*%e A (-'r*'y*'s+%i*'s*'x) , ' r) ; 

(d7) [r = - 1, r = 1] 

(c8) remvalue (fq) $ 

(c9) /* it should two identical real roots */ 
if length (ql) #2 then 'errorl$ 

(clO) if rhs (ql [1] ) #-rhs (ql [2] ) then 'error2$ 

(ell) r : abs (rhs (ql [ 1] ) ) ; 

(dll) 1 

(cl2) remvalue (ql) $ 

(cl3) /* the solution of the ODE at y>0 */ 
f xyp : ( ' ccl+' cc2 *' y) *%e A (-' 'r*abs {' s) * ' y-%i*' s * ' x) ; 

- abs (s) y - %i s x 

(dl3) (cc2 y + ccl) %e 

(cl4) /* the solution of the ODE at y<0 */ 

f xym: ( f cc3+'cc4* , y) *%e A ( ' ' r *abs ( ' s) *'y-%i*' s*'x) ; 

abs(s) y - %i s x 

( dl 4 ) (cc4 y + cc3) %e 

(cl5) /* calculate gsxxp */ 

gsxxp :dif f (fxyp , ' y, 2) ; 

2 - abs (s) y - %i s x 

(dl5) s (cc2 y + ccl) %e 

- abs (s) y - %i s x 

- 2 cc2 abs(s) %e 

(cl6) /* calculate gsyyp */ 
gsyyp :diff (fxyp, ' x, 2) ; 



- abs (s) y - %i s x 


2 

(dl6) - s (cc2 y + ccl) %e 

(cl7) /* calculate gtxyp */ 
gtxyp : dif f (fxyp, 'x, 1, 'y,l) ; 

- abs(s) y - %i s x 

(dl7) %i s abs(s) (cc2 y + ccl) %e 

- abs (s ) y - %i s x 

- %i cc2 s %e 

(cl8) /* calculate gsxxm */ 
gsxxm : dif f ( f xym, ' y , 2 ) ; 

2 abs(s) y - %i s x abs(s) y - %i s x 

(dl8) s (cc4 y + cc3) %e +2 cc4 abs(s) %e 

(cl9) /* calculate gsyym */ 

gsyym:diff (fxym, ' x, 2) ; 

2 abs(s) y - %i s x 

(dl9) - s (cc4 y + cc3) %e 

(c20) /* calculate gtxym */ 

gtxym:dif f (fxym, ' x, 1, 'y, 1) ; 

abs (s) y - %i s x 

(d20) - %i s abs(s) (cc4 y + cc3) %e 

abs (s) y - %i s x 

- %i cc4 s %e 

(c21) /* the shear tractions q(x)=0 so gtxy=0 for y=0 . 

According to this boundary condition cc2 can be expressed by ccl */ 

/* 

ui : gtxyp* %e A (%i*' s*' x) ; 
ui : rat subs t (0 , ' y, ui) ; 

11s : solve (ui, ' cc2 ) ; 
cc2 :rhs (11s [1] ) ; 

*/ 

/* Accirding to continuity conditions cc3 and cc4 could be expressed by ccl */ 

pp:diff(('ccl+'cc2*'y)*%e A (-'' r*abs ('s)* f y)-('cc3+'cc4*'y)*%e A (''r*abs('s)*'y) , ' 

abs(s) y abs(s) y 

(d21) - abs(s) (cc4 y + cc3) %e - cc4 %e 

- abs(s) y - abs(s) y 

- abs(s) (cc2 y + ccl) %e + cc2 %e 

(c22) ppp : ('ccl + ' cc2*'y) *%e A (-r*abs (s) *' y) - (' cc3+' cc4*' y) *%e A (r*abs (' s) *'y) ; 

- abs (s) y abs (s) y 

(d22) (cc2 y + ccl) %e - (cc4 y + cc3) %e 

(c23) ppp: rat subs t (0, 'y,ppp) ; 


(d23) 


ccl - cc3 



(c24) 11s : solve (ppp , ' cc3) ; 

(d24) [cc3 = cel] 

(c25) remvalue (ppp) $ 

(c26) cc3 :rhs (11s [1] ) ; 

(d26) cel 

(c27) ppt : pp ; 

abs (s) y abs (s) y 

(d27) - abs(s) (cc4 y + cc 3) %e - cc4 %e 


- abs (s ) y 

- abs(s) (cc2 y + ccl) %e + cc2 %e 


(c28) 

ppt :ratsubst (0, 'y,ppt) 

f 


(d28 ) 

(- cc3 - ccl) abs(s) - cc4 

+ cc2 

(c29) 

11s : solve (ppt , ' cc4 ) ; 



(d29) 

[cc4 = 

-- (- cc3 - ccl) abs(s) 

+ cc2] 

(c30) 

cc4 :rhs (11s [1] ) ; 



(d30) 

(- 

■ cc3 - ccl) abs (s) + 

cc2 

(c3 1 ) 

cc4 : ' ' cc4 ; 



(d31) 


cc2 - 2 ccl abs(s) 


(c32) 

remvalue (11s) $ 



(c33) 

remvalue (pp) $ 



(c34) 

remvalue (ppt) $ 



(c35) 

/* calculate fl (x) */ 



f bp : ( 

' ccl+cc2*' y) *%e A (- ,, r*abs( , s)*'y-%i*'s* , x) ; 


(d35) 

(cc2 y 

- abs(s) y 

■ + ccl) %e 

- %i s x 

(c36) 

fbm: (cc3+cc4*' y) *%e A ( ' 

'r*abs('s)*'y-%i*'s*' 

x) ; 

(d36) 

( (cc2 - 2 ccl 

abs(s) y - %i 

abs (s) ) y + ccl) %e 

(c37) 

eqql : integrate (fbp, ' x) 

/ 


(d37) 

%i (cc2 

- abs (s) y 

y + ccl) %e 

- %i s x 


s 


(c38) 

eqql : diff (eqql, ' y, 2 ) ; 



(d38 ) 

%i s (cc2 y + ccl) %e 

abs(s) y - %i s x 



2 %i cc2 abs(s) %e 


abs(s) y 


- abs (s) y - %i s x 



s 


(c39) eqq2 : dif f (f bp, ' x, 1) ; 


- abs (s) y - %i s x 


(d39) - %i s (cc2 y + ccl) %e 

(c40) f f 1 : ' all*eqql+' al2*eqq2 ; 

- abs (s ) y - %i s x 

(d40) all (%i s (cc2 y + ccl) %e 

abs (s) y - %i s x 


2 %i cc2 abs (s) %e 


-) 


- abs (s) y - %i s x 

- %i al2 s (cc2 y + ccl) %e 
(c41) ffl :ratsubst (0, 'y, ffl) ; 

2 - %i s x 

(2 %i all cc2 abs(s) + (%i al2 - %i all) ccl s ) %e 

(d41) 

s 

(c42) flxp:diff (ffl, ' x, 1) ; 

2 - %i s x 

(d42) %i (2 %i all cc2 abs(s) + (%i al2 - %i all) ccl s ) %e 
(c43) eqml : integrate (fbm, ' x) ; 

abs(s) y - %i s x 

%i ( (cc2 - 2 ccl abs(s)) y + ccl) %e 

(d43) 

s 

(c44) eqml : dif f (eqml, ' y, 2) ; 

abs (s) y - %i s x 

(d44) %i s ((cc2 - 2 ccl abs{s)) y + ccl) %e 

abs (s) y 

2 %i abs(s) (cc2 - 2 ccl abs(s)) %e 


(c45) eqm2 : dif f (fbm, ' x, 1) ; 

abs (s) y - %i s x 

(d45) - %i s ( (cc2 - 2 ccl abs(s)) y + ccl) %e 

(c46) fml : ' all*eqml+' al2*eqm2 ; 

abs(s) y - %i s x 

(d46) all (%i s ( (cc2 - 2 ccl abs(s)) y + ccl) %e 

abs (s) y - %i s x 

2 %i abs(s) (cc2 - 2 ccl abs(s)) %e 

+ ) 

s 


- %i s x 


abs (s) y - %i s x 



- %i al2 s ( (cc 2 - 2 ccl abs(s)) y + ccl) %e 
(c47) fml : ratsubst (0, ' y, fml) ; 

2 2 

(d47) - (- 2 %i all cc2 abs(s) + (%i al2 - %i all) ccl s + 4 %i all ccl s ) 

- %i s x 
%e /s 

(c48) f lxm: dif f (fml , ' x, 1) ; 

2 2 

(d48) %i (- 2 %i all cc2 abs(s) + (%i al2 - %i all) ccl s + 4 %i all ccl s ) 

- %i s x 
%e 

(c49) remvalue (eqql) $ 

(c50) remvalue (eqq2 ) $ 

(c51) remvalue (ff 2) $ 

(c52) remvalue (eqml) $ 

(c53) remvalue (eqm2) $ 

(c54) remvalue (fml) $ 

(c55) remvalue (ui) $ 

(c56) /* calculate fl (x) */ 
fix: (flxp-f lxm) / (2*%pi) ; 

2 - %i s x 

(d56) (%i (2 %i all cc2 abs(s) + (%i al2 - %i all) ccl s ) %e 

2 2 

- %i (-2 %i all cc2 abs(s) + (%i al2 - %i all) ccl s + 4 %i all ccl s ) 

- %i s x 

%e ) / (2 %pi) 

(c57) f lx : rat simp (fix) ; 

2 - %i s x 

(2 all cc2 abs(s) - 2 all ccl s ) %e 

(d 57) 

%pi 

(c58) mi : f lx*2*%pi*%e A (%i*s*x) ; 

2 

(d58) - 2 (2 all cc2 abs(s) - 2 all ccl s ) 

(c59) lx: solve (mi=f 1 (t)*%e A (%i*'s*'t) , ' cc2) ; 

%i s t 2 

%e fl(t) - 4 all ccl s 

(d59) ( cc2 ] 

4 all abs(s) 

(c60 ) /* 

lx : solve (mi=' integrate (f 1 (t ) *%e A (%i*' s*' t ) , t , -a, a) , ' cc2 ) ; 



*/ 

cc2 : rhs (lx [ 1] ) ; 

%i s t 2 

%e fl(t) - 4 all ccl s 

(cL60) 

4 all abs(s) 

(c61) /*ccl : rat subs t (f2 (t)*%e / '(%i*'s*'t) , ' integrate ( , %e /v (%i*'s*'t)*f2(t),t,-a,; 
f lx : ' integrate (f lx, s, mint , inf ) ; 

mi : ratsubst (2*all*cc2*abs ('s)-2*all*s A 2*cll,' integrate ( ('2*'all*'cc2*abs('s)-': 
*/ 

/* calculate f2 (x) */ 
eqql : integrate (fbp, ' y ) ; 

- abs(s) y - %i s x - abs(s) y - %i s x 

cc2 (abs(s) y + 1) %e ccl %e 

(d61) 

2 abs (s) 

s 

(c62) eqql : diff (eqql, ' x, 2) ; 

2 - abs (s) y - %i s x 

- abs(s) y - %i s x ccl s %e 

(d62) cc2 (abs(s) y + 1) %e + 

abs (s) 

(c63) eqq2 : diff (fbp , ' y, 1) ; 

- abs(s) y - %i s x - abs(s) y - %i s x 

(d63) cc2 %e - abs(s) (cc2 y + ccl) %e 

(c64) f f 2 : ' all*eqql+' al2*eqq2 ; 

- abs (s) y - %i s x 

(d64) all (cc2 (abs (s) y + 1) %e 
2 - abs (s) y - %i s x 

ccl s %e - abs(s) y - %i s x 

+ ) + al2 (cc2 %e 

abs(s) 

- abs (s) y - %i s x 
- abs(s) (cc2 y + ccl) %e ) 

(c65) ff2 :ratsubst (0, 'y, f f2) ; 

2 2 - %i s x 

( (- al2 - all) cc2 abs(s) + al2 ccl s - all ccl s ) %e 

(d65) 

abs (s) 

(c66) f 2xp : diff (ff2, ' x, 1) ; 

2 2 - %i s x 

%i s ( (- al2 - all) cc2 abs(s) + al2 ccl s - all ccl s ) %e 

(d66) 

abs (s) 


) , cc: 


*' ai: 



(c67 ) f 2xp : rat subs t ( ' ks*abs ( , s),'s / '3,f2xp); 


- -6 1 S X 

(d67 ) - ( (%i al2 + %i all) cc2 s + ccl (%i all ks - %i al2 ks)) %e 
(c68 ) f 2xp : rat sub st ( ' s*abs ( ' s ) , ' ks, f 2xp) ; 

(d68) - s (ccl (%i all abs(s) - %i al2 abs(s)) + (%i al2 + %i all) cc2) 

- %i s x 
%e 

(c69) eqml : integrate (fbm, ' y) ; 

abs (s) y - %i s x 

2 ccl abs(s) (abs(s) y - 1) %e 

(d69) 

2 

s 

abs(s) y - %i s x abs(s) y - %i s x 

cc2 (abs(s) y - 1) %e ccl %e 

+ + 

2 abs (s) 

s 

(c70) eqml : dif f (eqml, ' x, 2 ) ; 

abs(s) y - %i s x 

(d7Q) 2 ccl abs(s) (abs(s) y - 1) %e 

2 abs (s) y - %i s x 
abs(s) y - %i s x ccl s %e 

- cc2 (abs(s) y - 1) %e - 

abs (s) 

(c71) eqm2 : dif f (fbm, ' y, 1) ; 

abs (s ) y - %i s x 

(d71) abs(s) ((cc2 - 2 ccl abs(s)) y + ccl) %e 

abs (s) y - %i s x 

+ (cc2 - 2 ccl abs(s)) %e 

(c72) fm2 : ' all*eqml+' al2*eqm2 ; 

abs (s) y - %i s x 

(d72) al2 (abs(s) ( (cc2 - 2 ccl abs(s)) y + ccl) %e 

abs(s) y - %i s x 
+ (cc2 - 2 ccl abs(s)) %e ) 

abs ( s ) y - %i s x 

+ all (2 ccl abs(s) (abs(s) y - 1) %e 

2 abs(s) y-%isx 
abs(s) y - %i s x ccl s %e 

- cc2 (abs(s) y - 1) %e ) 

abs (s ) 

(c73) fm2 :ratsubst (0, 'y, fm2) ; 

(d73) 

2 2 - %i s x 

((- al2 - all) cc2 abs(s) + (al2 + 2 all) ccl s + all ccl s ) %e 



abs (s) 

(c74) f 2xm:diff (fm2, ' x, 1) ; 

(d74) 

2 2 - %i s x 

%i s ( (- al2 - all) cc2 abs(s) + (al2 + 2 all) ccl s + all ccl s ) %e 

abs (s) 

(c75) f 2xm: ratsubst ( ' ks*abs ('s),'s A 3,f2xm); 

- %i s x 

(d75) - ( (%i al2 + %i all) cc2 s + ccl (- %i al2 ks - 3 %i all ks)) %e 
(c76) f 2xm: ratsubst ( ' s*abs ( ' s ) , ' ks, f 2xm) ; 

(d76) - s (ccl (- %i al2 abs(s) - 3 %i all abs(s)) + (%i al2 + %i all) cc2) 

- %i s x 
%e 

(c77) remvalue (fxyp) $ 

(c78) remvalue (fxym) $ 

(c79) remvalue (r)$ 

(c80) remvalue (fbp) $ 

(c81) remvalue (fbm) $ 

(c82) remvalue (eqql) $ 

(c83) remvalue (eqq2 ) $ 

(c84) remvalue (ff 2) $ 

(c85) remvalue (eqml) $ 

(c86) remvalue (eqm2 ) $ 

(c87) remvalue (fm2) $ 

(c88) remvalue (ui) $ 

(c89) /* calculate f2(x) */ 
f 2x : ( f 2xp-f 2xm) / (2*%pi ) ; 

(d89) (s (ccl (- %i a!2 abs(s) - 3 %i all abs(s)) + (%i al2 + %i all) cc2) 

- %i s x 

%e - s (ccl (%i all abs(s) - %i al2 abs(s)) + (%i al2 + %i all) cc2) 

- %i s x 

%e ) / (2 %pi) 

(c90) f2x:ratsimp (f2x) ; 

- %i s x 

2 %i all ccl s abs(s) %e 


(d90) 


%pi 


(c91) f2x : ' integrate (f2x, s, mint, inf) ; 


inf 

/ 

[ - %i s x 

2 %i all ccl I s abs(s) %e ds 

] 

/ 

minf 


(d91) 


%pi 

(c92) /* solve ccl expressed by the auxiliary function */ 

mi : ratsubst ( ' s*abs ( ' s) , ' integrate {'s*abs('s)*%e A (-%i*s*x), s,minf , inf ) , f 2x) ; 

2 %i all ccl s abs(s) 

(d92) 


%pi 

(c93) mi:mi*2*%pi; 

(d93) - 4 %i all ccl s abs(s) 

(c94 ) lx : solve (mi=' integrate (f2 (t) *%e A <%i*' s* ' t ) , t , -a, a) , ' ccl) ; 


a 

/ 

[ %i s t 

%i I %e f 2 (t) dt 

] 

/ 

- a 

(d94) [ccl ] 

4 all s abs(s) 

(c95) ccl : rhs (lx [ 1] ) ; 


a 

/ 

[ %i s t 

%i I %e f 2 (t ) dt 

] 

/ 

- a 

(d95) 

4 all s abs (s) 

(c96) ccl : ratsubst (f2(t)*%e A (%i*'s*'t),'’ integrate ( , %e A (%i*'s*'t)*f2 (t) ,t,-a, 

%i s t 

%i %e f 2 (t ) 

(d96) 

4 all s abs (s) 

(c97) /* calculate cc2 cc3 and cc4 expressed by the auxiliary function */ 
cc2 : ' ' ccl; 


(d97) 


%i s t 

%i s t %i s %e f 2 (t ) 

%e fl(t) 

abs (s) 


4 all abs(s) 


) , ccl) 



(c98) cc2 : xthru (cc2 ) ; 


(d98) 

(c99) 

(d99) 

(clOO) 

(dlOO ) 

(clOl ) 

(dlOl) 

(cl02) 

(dl02) 
(cl03 ) 

(dl03) 
(cl04 ) 

(dl04 ) 
(cl05) 

(dl05) 

(cl06) 


%i s t %i s t 

abs(s) %e fl(t) - %i s %e f2(t) 


4 all s 


cc2 : expand (cc2 ) ; 


%i s t 


;i s t 


%i %e f2(t> abs(s) %e fl(t) 


4 all s 


4 all s 


cc21 : last (cc2 ) ; 


%i s t 

abs (s) %e f 1 (t) 


4 all s 


cc2 : first (cc2 ) ; 


%i s t 
4 all s 


%i %e f 2 (t) 


cc21 : ratsubst (kb,s A 2,cc21) ; 


%i s t 

abs (s) %e f 1 (t) 

4 all kb 

cc21 : ratsubst (1, abs (s) , cc21) ; 

%i s t 

%e fl(t) 


4 all kb 

cc21 :ratsubst (abs (s) , kb, cc21) ; 

%i s t 

%e fl(t) 

4 all abs (s) 


cc2 : cc2+cc21; 


%i s t 


%i s t 


%i %e f 2 (t) %e fl (t) 


4 all s 


cc3 : ' ' cc3; 


%i s t 


si %e 


4 all abs(s) 


f 2 (t ) 


(did 6) 


4 all s abs(s) 



(cl07 ) cc4 : ' ' cc4 ; 


%i s t %i s t 

%i %e f 2 (t) %e fl (t) 

( dl 07) 

4 all s 4 all abs (s) 

(cl08) /* Get the answer of gsxxp expressed by the auxiliary function */ 
mwp : ' ' gsxxp; 


%i s t %i s t %i s t 

2 %i %e f 2 (t ) %e fl(t) %i %e f2(t) 

( dl 08) s (( ) y + ) 

4 all s 4 all abs(s) 4 all s abs(s) 

%i s t %i s t 

- abs(s) y - %i s x %i %e f2 (t) %e fl(t) 

%e - 2 abs(s) ( ) 

4 all s 4 all abs (s) 


- abs(s) y - %i s x 
%e 

(cl09) mwp : expand (mwp) ; 

- abs(s) y-%isx+%ist 

%i s f 2 (t ) y %e 

(dl09) 

4 all 

2 - abs(s) y-%isx+%ist 

s fl(t) y %e 


4 all abs(s) 

- abs(s) y-%isx+%ist 

%i abs (s ) f 2 (t ) %e 


2 all s 


- abs(s) y-%isx+%ist 

%i s f2 (t) %e 

+ 

4 all abs(s) 

- abs(s) y - %i s x + %i s t 

fl(t) %e 

+ 

2 all 

(cllO) /* calculate |s| when s>0 */ 
assume (s>0) ; 


(dllO ) 
(clll) 

(dill) 


[s > 0] 


mip : ' 'mwp; 

-sy-%isx+%ist 

%i s f 2 (t ) y %e 

4 all 


sy-%isx+%ist 


s y 


%i s x + %i s t 



s f 1 (t ) y %e 


%i f 2 (t ) %e 


4 all 


4 all 


fl(t) %e 


-sy-%isx+%ist 


2 all 


(cll2) mip : expand (mip) ; 


(dll2 ) 


%i s f 2 (t ) y %e 


-sy-%isx+%ist 


4 all 


-sy-%isx+%ist -sy-%isx+%ist 

s fl(t) y %e %i f2(t) %e 


4 all 


4 all 


fl(t) %e 


-sy-%isx+%ist 


2 all 


(cll3) forget (s>0); 
(dll3) 


[s > 0] 


(cll4) /* calculate |s| when s<0 */ 
assume (s<0) ; 


(dll4 ) 

(cll5) mim:''mwp; 


[s < 0] 


sy-%isx+%ist sy-%isx+%ist 

%i s f2(t) y %e s fl(t) y %e 

(dll5 ) + 

4 all 4 all 


sy-%isx+%ist sy-%isx+%ist 

%i f 2 (t) %e fl(t) %e 


4 all 


2 all 


(cll6) mim: ratsubst (-s, s,mim) ; 

%i s x %i s x 

(dll6) - ( (%i s f 2 (t ) + s fl(t)) %e y + (- %i f2(t) - 2 fl(t)) %e ) 


- s y - %i s t 


/ (4 all) 


[cll7) mim: expand (mim) ; 


(dll7) - 


)i s f 2 (t ) y %e 


-sy+%isx-%ist 


4 all 


-sy+%isx-%ist -sy+%isx-%ist 

s fl(t) y %e %i f 2 (t ) %e 

+ 



4 all 


4 all 


-sy+%isx-%ist 

fl(t) %e 

+ 

2 all 

(cll8) forget (s<0); 

(dll8) [s < 0] 

(cll9) /* combine the two terms */ 
mi : mip+mim; 

-sy+%isx-%ist 


(dll9 ) - 


%i s f 2 (t ) y %e 


4 all 


-sy+%isx-%ist -sy+%isx-%ist 

s fl(t) y %e %i f2(t) %e 

+ 

4 all 4 all 

-sy+%isx-%ist -sy-%isx+%ist 

fl (t) %e %i s f 2 (t) y %e 

+ + 

2 all 4 all 

-sy-%isx+%ist -sy-%isx+%ist 

s fl(t) y %e %i f2(t) %e 


fl(t) %e 


4 all 

-sy-%isx+%ist 


4 all 


2 all 

(cl 20 ) mi : subst (wp, %e A (-s*y-%i*s*x+%i*s*t ) , mi) ; 

-sy+%isx-%ist 


(dl 20) - 


%i s f 2 (t ) y %e 


4 all 

-sy+%isx-%ist 

s f 1 (t ) y %e 

4 all 

-sy+%isx-%ist 


;i f2(t) %e 


-sy+%isx-%ist 


4 all 


fl(t) %e 


%i s f2(t) wp y s fl(t) wp y 


2 all 

%i f2(t) wp fl(t) wp 

+ 

4 all 2 all 


4 all 


4 all 


(cl21) mi : subst (wm, %e A (-s*y+%i*s*x-%i*s*t ) ,mi) ; 

%i s f2(t) wp y s fl(t) wp y %i s f2(t) wm y s fl(t) wm y 


(d!21) 


4 all 


4 all 


4 . all 


4 all 




{ cl22 ) 
(dl22 ) 

(cl23) 

(dl23) 

(cl24 ) 
(dl24) 

(cl25) 
(dl25 ) 

(cl26) 

(dl26) 

(cl27 ) 
(dl27 ) 

(cl28 ) 
(dl28 ) 
(cl29) 
(dl29 ) 

(cl30 ) 


%i f 2 (t ) wp fl(t) wp %i f2(t) wm fl(t) wm 

— + + + 

4 all 2 all 4 all 2 all 

mi :ratsubst (ywp, s*wp,mi) ; 

- (y (- %i f2(t) ywp + fl(t) ywp + (%i s f2(t) + s fl(t)) wm) 

+ (%i f2(t) - 2 fl(t)) wp + (- %i f 2 (t ) - 2 fl(t)) wm)/(4 all) 
mi : ratsubst (ywm, s*wm,mi) ; 

((%i f2(t) - fl(t)) y ywp + y (- %i f2(t) ywm - fl(t) ywm) 

+ (2 fl(t) - %i f 2 (t ) ) wp + (%i f 2 (t ) + 2 fl(t)) wm) / ( 4 all) 

ywp: (1/ ( (-y+%i* (t-x) ) A 2) ) ; 

1 

2 

(%i (t-x) - y) 

mi : mi ; 

( (%i f2(t) - fl(t)) y ywp + y (- %i f2(t) ywm - fl(t) ywm) 

+ (2 fl(t) - %i f 2 (t ) ) wp + (%i f 2 (t ) + 2 fl(t)) wm)/(4 all) 

ywm: (1/ ( (-y-%i* (t-x) ) A 2) ) ; 

1 


(- y - %i (t-x)) 


(%i f 2 (t) - fl (t) ) y %i f 2 (t) 

( + ( 

2 2 
(%i (t - x) - y) (- y - %i (t - x) ) 


fl (t) 


■) y 


(- y - %i (t - x) ) 

+ (2 fl(t) - %i f 2 (t ) ) wp + (%i f 2 (t ) + 2 fl(t)) wm)/(4 all) 
wp : (-1/ (-y+%i* (t-x) ) ) ; 

1 




i (t - x) - y 


mx : mx ; 

(%i f 2 (t) - fl(t)) y 


(- 


%i f 2 (t) fl(t) 

___ + ( ) y 

2 2 2 
(- y - %i (t - x) ) (- y - %i (t - x) ) 


(%i (t-x) - y) 

+ (2 fl(t) - %i f 2 (t ) ) wp + (%i f 2 (t ) + 2 fl(t)) wm)/(4 all) 
wm: (-1/ (-y-%i* (t-x) ) ) ; 

1 


(dl30 ) 


y 


%i (t - x) 



( cl31 ) 


mi : ' ' mi ; 


(%i f 2 (t ) - fl(t)) y %i f 2 (t ) fl(t) 

(d!31 ) ( + ( 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x) ) (- y - %i (t - x) ) 


) y 


2 fl(t> - %i f 2 (t) %i f 2 (t ) + 2 fl{t) 


(cl 32) mi :mi/ (2*%pi) ; 


%i (t - x) - y 


%i f 2 (t) 


- y - %i (t - x) 


fl (t) 


■)/(4 all) 


(%i f2(t) - fl(t)) y 

(dl32 ) ( + ( ; 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x) ) (- y - %i (t - x) ) 


2 fl(t) - %i f 2 (t ) %i f 2 (t ) + 2 fl(t) 


%i (t - x) - y 


- y - %i (t - x) 


-) / (8 %pi all) 


(cl33) mi : subst (z , (t-x) ,mi) ; 

%i f 2 (t) 


(dl33) (y (- 


fl (t) 2 fl (t) - %i f 2 (t) 

) 

2 2 %i z - y 

(- %i z - y) (- %i z - y) 

(%i f 2 (t ) - fl(t)) y %i f 2 (t ) + 2 fl(t) 

+ ) / (8 % p i all) 

2 - %i z - y 

(%i z - y) 

(cl34) mi : distrib (mi) ; 

/usr/local/macsyma_30 9/mrg/scs . o being loaded. 

%i f 2 (t) fl (t) 

y ( ) 

2 2 
(- %i z - y) (- %i z - y) 


(dl34 ) 


8 %pi all 


2 fl (t) - %i f 2 (t) 

8 %pi all (%i z - y) 


(%i f 2 (t ) - fl(t)) y %i f2(t) + 2 fl(t) 


2 8 %pi all (- %i z - y) 

8 %pi all (%i z - y) 

(cl35) ki : first (mi) ; 


(dl35) 


%i f 2 (t) fl (t) 

y ( ) 

2 2 
(- %i z - y) (- %i z - y) 


8 %pi all 


(cl36) kil : rest (mi, 2) ; 


(dl36 ) 


(%i f 2 (t) - fl(t)) y %i f 2 (t ) + 2 fl(t) 


2 8 %pi all (- %i z - y) 

8 %pi all (%i z - y) 



(cl37) ki : ki+f irst (kil) ; 


%i f 2 {t ) fl(t) 

y { ) 

2 2 

(- %i z - y) (- %i z - y) (%i f2(t) - fl(t)) y 

(dl37) + 

8 %pi all 2 

8 %pi all (%i z - y) 

(cl38) ki : xthru (ki) ; 

2 2 
(- %i f 2 (t ) - fl(t)) y (%i z - y) + (%i f2(t) - fl(t)) y (- %i z - y) 

(dl38 ) 

2 2 

8 %pi all (- %i z - y) (%i z - y) 

(cl39) ki : expand (ki) ; 

2 

2 fl(t) y z 

(dl39 ) 

4 2 2 4 

8 %pi all z +16 %pi all y z + 8 %pi all y 

2 

4 f2(t) y z 


4 2 2 4 

8 %pi all z +16 %pi all y z + 8 %pi all y 

3 

2 fl(t) y 


4 2 2 4 

8 %pi all z +16 %pi all y z + 8 %pi all y 

(cl40) ki : factor (ki); 


(dl40 ) 


2 2 
y (fl(t) z - 2 f 2 (t ) y z - fl(t) y ) 

2 2 2 

4 %pi all (z + y ) 


(cl41) k j : rest (mi, 1) ; 

2 fl (t) - %i f 2 (t) 


(%i f 2 (t) - fl(t)) y 


(dl41) 


8 %pi all (%i z - y) 


8 %pi all (%i z - y) 


%i f 2 (t ) +2 fl (t) 

8 %pi all (- %i z - y) 


( cl42 ) 
(dl42 ) 
(cl 43) 
(dl43) 


k j : first (kj ) ; 


2 fl(t) - %i f 2 (t ) 

8 %pi all (%i z - y) 


k j : k j+last (mi) ; 

2 fl(t) - %i f 2 (t) %i f 2 (t ) + 2 fl(t) 

8 %pi all (%i z - y) 8 %pi all (- %i z - y) 



(cl44) k j : xthru (k j ) ; 


(d.14 4 ) 
(Cl45) 
(dl45) 

(cl46) 

(dl46) 

(Cl47 ) 

(dl47) 

(cl48 ) 

(dl48) 

(cl49) 
gsxxp : 

(dl49) 

(C150) 


(- %i f 2 (t ) - 2 fl(t)) ( %i z - y) + (%i f2(t) - 2 fl(t)) (- %i z - y) 


8 %pi all (- %i z - y) (%i z - y) 

k j : expand (k j ) ; 

2 f 2 (t ) z 4 fl(t) y 


2 2 2 2 
8 %pi all z +8 %pi all y 8 %pi all z +8 %pi all y 

k j : factor (k j ) ; 

f 2 (t) z + 2 fl (t) y 


2 2 

4 %pi all (z + y ) 

mi : ki+k j ; 

2 2 

y (fl(t) z - 2 f2(t) y z - fl(t) y ) f2(t) z + 2 fl(t) y 


2 2 2 2 2 
4 %pi all (z + y ) 4 %pi all (z + y ) 

mi : subst ( (t-x) , z,mi) ; 

2 2 
y (- fl(t) y - 2 f 2 (t ) (t-x) y + fl(t) (t-x) ) 


2 2 2 
4 %pi all (y + (t - x) ) 

2 fl(t) y + f2 (t) (t - x) 


2 2 
4 %pi all (y + (t - x) ) 

/* Get the answer of gsxxp expressed by the auxiliary function */ 
integrate (mi, t , -a, a) ; 
a 

/ 2 2 
[ y (- fl(t) y - 2 f 2 (t ) (t-x) y + fl(t) (t-x) ) 

I ( 

] 2 2 2 
/ 4 %pi all (y + (t - x) ) 

- a 

2 fl (t) y + f 2 (t) (t - x) 

+ ) dt 

2 2 
4 %pi all (y + (t - x) ) 

ki:ki*4*%pi*all; 

2 2 
y (fl(t) z - 2 f2(t) y z - fl(t) y ) 


2 2 2 
(z + y ) 


(dl50) 



(cl51) kj:kj*4*%pi*all; 


(dl51 ) 


f 2 (t ) z + 2 fl(t) y 

2 2 
z + y 


(cl52) mi:ki+kj; 

2 2 

y (fl(t) z - 2 f2(t) y z - fl(t) y ) f2(t) z + 2 fl(t) y 

(dl52 ) + 

2 2 2 2 2 
(z + y ) z + y 

(cl53) gsxxpf : subst ( (t-x) , z,mi) ; 

2 2 
y (- fl(t) y - 2 f 2 (t ) (t-x) y + fl(t) (t-x) ) 


(d!53 ) 


2 2 2 
(y + (t - x) ) 


2 fl(t) y + f2 (t) (t - x) 

+ 

2 2 
y + (t - x) 

(cl54) /* calculate the gsyyp and make it simplified */ 
mwp : ' ' gsyyp; 

%i s t %i s t %i s t 

2 %i %e f2(t) %e fl(t) %i %e f2(t) 

(dl54 ) - s (( ) y + — ) 

4 all s 4 all abs(s) 4 all s abs(s) 

- abs (s ) y - %i s x 


(cl55) mwp : expand (mwp) ; 


(dl55) - 


%i s f2 (t) y %e 


- abs(s) y-%isx+%ist 


4 all 


2 - abs(s) y-%isx+%ist 

s fl(t) y %e 


%i s f 2 (t ) %e 


4 all abs(s) 

• abs(s) y-%isx+%ist 


4 all abs(s) 

(cl56) /* calculate |s| when s>0 */ 


assume (s>0) ; 
(dl56) 


[s > 0] 


(cl 57) mip ; ' ' mwp ; 


-sy-%isx+%ist 



%i s f 2 ( t ) y %e 

(dl57 ) 

4 all 

-sy-%isx+%ist -sy-%isx+%ist 

s fl(t) y %e %i f2(t) %e 

+ 

4 all 4 all 

(cl58) mip:expand(mip); 

-sy-%isx+%ist 

%i s f 2 (t ) y %e 

(dl58 ) 

4 all 

-sy-%isx+%ist -sy-%isx+%ist 

s fl(t) y %e %i f2 (t) %e 

+ 

4 all 4 all 


(cl59) forget (s>0); 

(dl59) [s > 0] 

(cl60) /* calculate |s| when s<0 */ 
assume (s<0) ; 


(dl60) 

(cl61) 

(dl61) 


[s < 0] 

mim : ' ' mwp ; 

sy-%isx+%ist 

%i s f2 (t) y %e 


4 all 

sy-%isx+%ist sy-%isx+%ist 

s fl(t) y %e %i f2(t) %e 


4 all 4 all 


(cl62) mim: ratsubst (-s, s,mim) ; 

%isx %isx - s y - %i s t 

( (%i s f2(t) + s fl(t)) %e y + %i f2(t) %e ) %e 

<dl62) 

4 all 


(cl63) mim: expand (mim) ; 


(dl63) 


-sy+%isx-%ist 

%i s f2 (t) y %e 


4 all 


+ 


s f 1 (t) y %e 


sy+%isx-%ist 

+ 

4 all 


- s y + %i 

%i f 2 (t) %e 


4 all 


x - %i s t 


(cl64) forget (s<0); 


(dl64 ) 


[s < 0] 



(cl65) /* combine the two terms */ 
mi :mip+mim; 


(d!65) 


-sy + %isx-%ist 

%i s f 2 (t ) y %e 


4 all 


s fl(t) y %e 


-sy+%isx-%ist 


%i f 2 (t) %e 


- s y + 


% -i 


1 s x 


- %-i 


%i s t 


4 all 


4 all 


-sy-%isx+%ist 

%i s f 2 (t) y %e 


4 all 


-sy-%isx+%ist 

s fl(t) y %e 


4 all 


-sy-%isx+%ist 

%i f 2 (t ) %e 


4 all 

(cl 66) mi : subst (wwp, %e A (-s*y-%i*s*x+%i*s*t) ,mi) ; 

-sy+%isx-%ist 

%i s f2(t) y %e 

(dl66) 

4 all 

-sy+%isx-%ist -sy+%isx-%ist 

s fl(t) y %e %i f2(t) %e 

+ + 

4 all 4 all 

%i s f2(t) wwp y s fl(t) wwp y %i f2(t) wwp 


4 all 4 all 4 all 

(cl 67) mi : subst (wwm, %e A (-s*y+%i*s*x-%i*s*t) ,mi) ; 

%i s f2 (t) wwp y s fl(t) wwp y %i s f2(t) wwm y s fl (t) wwm y 

(dl67) + + + 

4 all 4 all 4 all 4 all 

%i f2 (t) wwp %i f2(t) wwm 


4 all 4 all 

(cl68) mi rratsubst (ywwp, s*wwp,mi) ; 

(dl68) (y (- %i f 2 (t) ywwp + fl(t) ywwp + (%i s f2(t) + s fl(t)) wwm) 

- %i f2(t) wwp + %i f2(t) wwm) / (4 all) 

(cl69) mi :ratsubst (ywwm, s*wwm,mi) ; 

(dl69) - ( (%i f2 (t) - fl(t)) y ywwp + y (- %i f2(t) ywwm - fl(t) ywwm) 

+ %i f2(t) wwp - %i f2(t) wwm)/ (4 all) 

(cl 70 ) ywwp : (1/ ( (-y+%i* (t-x) ) A 2) ) ; 

1 

(dl 70) 



2 


(%i (t - x) - y) 

( c 1 7 1 ) mi : mi ; 

(dl71) - ( (%i f2(t) - fl(t)) y ywwp + y (- %i f2(t) ywwm - fl(t) ywwm) 

+ %i f 2 (t ) wwp - %i f 2 (t ) wwm) / (4 all) 

(C172) ywwm: (1/ ( (-y-%i* (t-x) ) A 2) ) ; 


(d.172 ) 


(cl 73) mi : ' 'mi; 


1 


2 

(- y - %i (t - x) ) 


(%i f 2 (t ) - fl(t)) y %i f 2 (t ) fl(t) 

(dl73) - ( + ( 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x) ) (- y - %i (t - x) ) 


y + %i f2(t) wwp - %i f2 (t) wwm)/ (4 all) 


(cl74) wwp : (-1/ (-y+%i* (t-x) ) ) ; 
(dl74 ) 

(cl 75) mi : mi ; 


%i (t-x) - y 


%i f 2 (t ) 


fl (t) 


(%i f 2 (t ) - fl(t)) y 

(dl75 ) - ( + ( ] 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x)) (- y - %i (t-x)) 


y + %i f2(t) wwp - %i f2(t) wwm)/ (4 all) 


(cl76) wwm: (-1/ (-y-%i* (t-x) ) ) ; 
(dl76) 

(cl77) mi : ' 'mi; 


- y - %i (t-x) 


%i f 2 (t) 


fl (t) 


(%i f 2 (t ) - fl(t)) y 

(dl77 ) - ( + ( 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x) ) {- y - %i (t - x) ) 


%i f 2 (t) %i f 2 (t) 

y + ) / (4 a n) 

%i (t - x) - y - y - %i (t - x) 


(cl78) mi :mi/ (2*%pi) ; 


%i f 2 (t ) 


fl (t) 


(%i f 2 (t) - fl(t)) y 

(dl78) - ( + ( ) 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x) ) (- y - %i (t - x) ) 

%i f 2 (t) %i f 2 (t) 

y + ) / (8 %pi all) 



(C179) 

(dl79) 

(cl80) 

(dl80) 

( cl81 ) 

(dl81) 
(cl82 ) 
(dl82 ) 

(cl83) 

(dl83) 

(cl84 ) 
(dl84) 

(cl85 ) 


%i (t - x) - y - y - %i (t - x) 


mi : subst (z, (t-x) ,mi) ; 

%i f 2 (t) 

2 


fl(t) 


- (y (- 


(- %i z - y) (- %i z - y) 


%i f 2 (t ) (%i f2(t) - fl(t)) y 

+ 

%i z - y 2 

(%i z - y) 

%i f 2 (t) 

+ ) / (8 %pi all) 

- %i z - y 


mi :distrib (mi) ; 

%i f 2 (t) 


fl(t) 


y (- 


(- %i z - y) (- %i z - y) %i f2 (t) 

+ 

8 %pi all 8 %pi all (%i z - y) 

(%i f 2 (t) - fl (t) ) y %i f 2 (t) 

2 8 %pi all (- %i z - y) 


8 %pi all (%i z - y) 


ki : first (mi) ; 


%i f 2 (t) fl (t) 

y ( ) 

2 2 
(- %i z - y) (- %i z - y) 

8 %pi all 

kil : rest (mi, 2) ; 

(%i f 2 (t) - fl (t) ) y %i f 2 (t) 

2 8 %pi all (- %i z - y) 

8 %pi all (%i z - y) 
ki : ki+first (kil) ; 

%i f 2 (t) fl (t) 


y ( ) 

2 2 
(- %i z - y) (- %i z - y) 


(%i f 2 (t) - fl (t) ) y 


8 %pi all 


8 %pi all (%i z - y) 


ki : xthru (ki) ; 


2 2 
(- %i f2(t) - fl(t)) y (%i z - y) - (%i f2(t) - fl(t)) y (- %i z - y) 


2 2 

8 %pi all (- %i z - y) (%i z - y) 

ki : expand (ki) ; 


2 



2 f 1 (t ) y z 

(dl85) 

4 2 2 4 

8 %pi all z +16 %pi all y z +8 %pi all y 

2 

4 f 2 (t ) y z 

+ 

4 2 2 4 

8 %pi all z +16 %pi all y z +8 %pi all y 

3 

2 fl(t) y 

+ 

4 2 2 4 

8 %pi all z +16 %pi all y z + 8 %pi all y 

(cl86) ki : factor (ki) ; 

2 2 
y (fl(t) z - 2 f 2 (t) y z - fl(t) y ) 

(dl 86) 

2 2 2 

4 %pi all (z + y ) 

(cl87) k j : rest (mi, 1) ; 

%i f 2 (t ) (%i f 2 (t) - fl(t)) y %i 

(dl87) 

8 %pi all (%i z - y) 28 %pi all 

8 %pi all (%i z - y) 

(cl88) k j : first (kj ) ; 

%i f 2 (t) 

(dl88) 

8 %pi all (%i z - y) 

(cl89) k j : k j+last (mi) ; 

%i f 2 (t ) %i f2(t) 

(d 189 ) 

8 %pi all (%i z - y) 8 %pi all (- %i z - y) 

(cl90) k j : xthru (k j ) ; 

%i f2(t) (- %i z - y) - %i f2(t) (%i z - y) 

(di90 ) 

8 %pi all (- %i z - y) (%i z - y) 

(cl91) k j : expand (kj ) ; 

f2(t) z 

(di9i) 

2 2 
4 %pi all z +4 %pi all y 

(cl92) k j : factor (kj ) ; 

f 2 (t ) z 

(di 92 ) 

2 2 

4 %pi all (z + y ) 


f 2 (t) 

(- %i z - y) 


(cl93) mi:ki+kj; 



f 2 (t ) z 


(dl93) 


2 2 

4 %pi all (z + y ) 


(cl94) mi : subst ( (t-x) , z,mi) ; 


2 2 
y (fl(t) z - 2 f 2 (t ) y z - fl(t) y ) 


2 2 2 

4 %pi all (z + y ) 


(dl94 ) 


f 2 (t ) (t - x) 


2 2 
4 %pi all (y + (t - x) ) 


2 

y (- fl(t> y - 2 f 2 (t) (t-x) y + fl(t) (t 


2 2 2 
4 %pi all (y + (t - x) ) 

(cl95) /* Get the answer of gsyyp expressed by the auxiliary function */ 
gsyyp : ' integrate (mi, t , -a, a) ; 


2 

x) ) 


/ 

[ f 2 (t) (t - x) 

(dl95) I ( 

] 2 2 
/ 4 %pi all (y + (t - x) ) 

- a 


2 

y (- fl(t) y - 2 f2(t) (t-x) y + fl(t) 


2 2 2 
4 %pi all (y + (t - x) ) 


2 

(t-x) ) 
) dt 


(cl96) ki : ki*4*%pi*all; 


2 2 
y (fl(t) z - 2 f 2 (t ) y z - fl(t) y ) 

(dl96) 

2 2 2 

(z + y ) 

(cl97) kj :kj*4*%pi*all; 

f2(t) z 

(dl97 ) 

2 2 

z + y 

(cl98) mi:ki+kj; 

2 2 
f2(t) z y (fl(t) z - 2 f2(t) y z - fl(t) y ) 

(dl98 ) 

2 2 2 2 2 
z + y (z + y ) 

( cl 9 9 ) gsyypf :subst ( (t-x) , z,mi) ; 


f 2 (t ) (t - x) 


2 

y (- fl(t) y - 2 f 2 (t ) (t-x) y + fl(t) (t 


2 2 


2 

x) ) 


(dl99) 


2 


2 


2 



y + (t - x) (y + (t - x) ) 
(c200) /* calculate the gtxyp and make it simplified */ 
mwp : ' ' gtxyp; 


%i s t %i s t %i 

%i %e f 2 (t) %e fl (t) %i %e 

(d200) %i s abs(s) (( ) y + 

4 all s 4 all abs(s) 4 all 

%i s t %i s t 

- abs(s) y - %i s x %i %e f2(t) %e fl(t) 

%e - %i s ( ) 

4 all s 4 all abs (s) 


-abs(s) y-%isx 
%e 

(c201) mwp : expand (mwp) ; 

- abs(s) y-%isx+%ist 

abs(s) f2(t) y %e 

(d201 ) 

4 all 

- abs(s) y-%isx+%ist 

%i s fl(t) y %e 


4 all 

- abs(s) y-%isx+%ist 

%i s f 1 (t) %e 

+ 

4 all abs(s) 

(c202) /* calculate |s| when s>0 */ 
assume ( s>0 ) ; 

(d202 ) [s > 0] 

(c203) mip:''mwp; 

-sy-%isx+%ist 

s f2(t) y %e 

(d203) 

4 all 

-sy-%isx+%ist - s y - %i s 

%i s fl(t) y %e %i fl (t) %e 


4 all 4 all 


(c204) mip : expand (mip) ; 

-sy-%isx+%ist 


(d204 ) - 


s f 2 (t) y %e 


4 all 

-sy-%isx+%ist 


si s f 1 (t ) y %e 


%i f 1 (t) %e 


- s y - 


s t 

f 2 (t) 

) 

abs (s ) 


x + %i s t 


x + %i s t 


4 all 


+ 


4 all 



(c205) forget (s>0); 

(d205) [s > 0] 

(c206) /* calculate |s| when s<0 */ 
assume (s<0) ; 

(d206) [s < 0] 

(c207) mim: ' ' mwp; 

sy-%isx+%ist sy-%isx+%ist 

s f2(t) y %e %i s fl(t) y %e 

(d207 ) 

4 all 4 all 

sy-%isx+%ist 

%i fl(t) %e 

4 all 

(c208) mim: ratsubst (-s, s, mim) ; 

( d2 0 8 ) 

%isx %isx - s y - %i s t 

( (s f2(t) - %i s fl(t)) %e y + %i fl(t) %e ) %e 

4 all 


(c209) mim: expand (mim) ; 


(d2 09 ) 


-sy+%isx-%ist 

s f 2 (t) y %e 


4 all 


-sy+%isx-%ist 

%i s f 1 (t ) y %e 

+ 

4 all 


-sy+%isx-%ist 

%i fl(t) %e 


4 all 


(c210) forget (s<0); 

(d2 10 ) [s < 0] 

(c211) /* combine the two terms */ 


mi :mip+mim; 


(d211) 


-sy+%isx-%ist 

s f 2 (t ) y %e 


4 all 


%i s fl(t) y %e 


-sy+%isx-%ist 


%i fl (t) %e 


-sy+%isx-%ist 


s f2(t) y %e 


4 all 

-sy-%isx+%ist 


%i s fl(t) y %e 


4 all 

sy-%isx+%ist 


4 all 


4 all 



-sy-%isx+%ist 

%i fl(t) %e 

+ 

4 all 

(c212) mi : subst (twp, %e A (-s*y-%i*s*x+%i*s*t ) , mi) ; 

-sy+%isx-%ist 

s f2 (t ) y %e 

(d212) 

4 all 


-sy+%isx-%ist -sy+%isx-%ist 

%i s fl(t) y %e %i fl(t) %e 

+ 

4 all 4 all 

s f2(t) twp y %i s fl(t) twp y %i fl(t) twp 

+ 

4 all 4 all 4 all 

(c213) mi : subst (twm, %e A (-s*y+%i*s*x-%i*s*t) ,mi) ; 


s f2(t) twp y %i s fl(t) twp y s f2(t) twm y %i s fl (t) twm y 

(d213) + 

4 all 4 all 4 all 4 all 


%i fl(t) twp %i f 1 (t) twm 
4 all 4 all 


(c214 ) mi : rat subst (twwp, s*twp,mi) ; 

(d2l4) - ((f2(t) twwp + %i fl(t) twwp + (s f2(t) - %i s fl(t)) twm) y 

- %i fl(t) twp + %i fl(t) twm)/ (4 all) 

(c215) mi :ratsubst (twwm, s*twm,mi) ; 

(d215) - ( ( (f2 (t) + %i fl(t)) twwp + f2(t) twwm - %i fl(t) twwm) y 

- %i fl(t) twp + %i fl(t) twm)/ (4 all) 

(c216) twwp: (1/ ( {-y+%i* (t-x) ) A 2) ) ; 


(d216) 


1 


2 

(%i (t-x) - y) 


(c217) mi:mi; 

(d217) - (((f2(t) + %i fl(t)) twwp + f2 (t) twwm - %i fl(t) twwm) y 

- %i fl(t) twp + %i fl(t) twm)/ (4 all) 


(c2l8) twwm: (1/ ( (-y-%i* (t-x) ) A 2) ) ; 


(d2l8) 


1 


2 

(- y - %i (t - x) } 


(c2l9) mi : ' 'mi; 


f 2 (t) + %i fl (t) 


f 2 (t ) 


%i fl (t) 



(d219) - (( 


+ 


2 2 2 
(%i (t - x) - y) (- y - %i (t - x) ) (- y - %i (t - x)) 


■) Y 


- %i fl(t) twp + %i fl(t) twm) / ( 4 all) 


(c220) twp: (-1/ (-y+%i* (t-x) ) ) ; 
(d220) 

(c22l) mi: mi; 


%i (t-x) - y 


f 2 (t ) 


%i fl(t) 


f 2 (t) + %i fl (t) 

(d22l) - (( + 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x)) (- y - %i (t-x)) 


) y 


- %i fl(t) twp + %i fl(t) twm)/ (4 all) 


(c222) twm: (-1/ (-y-%i* (t-x) ) ) ; 
(d222 ) 

(c223) mi : ' 'mi; 


- y - %i (t-x) 


f 2 (t ) 


%i fl (t) 


f 2 (t) + %i fl(t) 

(d223) - (( + 

2 2 2 
(%i (t - x) - y) (- y - %i (t - x)) (- y - %i (t-x)) 


) y 


%i fl (t) %i fl (t) 

+ — ) / (4 all) 

%i (t - x) - y - y - %i (t - x) 


(c224) mi :mi/ (2*%pi) ; 

f 2 (t) + %i fl (t) 

( d2 2 4 ) - (( + 

2 


f 2 (t) 


%i fl(t) 


•) y 


(%i (t - x) - y) (- y - %i (t - x)) (- y - %i (t-x)) 


%i fl (t) %i fl (t) 

+ ) / ( 8 % p i all) 

%i (t - x) - y - y - %i (t - x) 


(c225) mi : subst (z , (t-x) ,mi) ; 

(d225) 

f 2 (t) + %i fl (t) f 2 (t) 

y ( + 

2 2 


%i fl (t) %i fl (t) %i fl (t) 

■) + 


2 %i z - y - %i z - y 


( %i z - y) 


(- %i z - y) (- %i z - y) 


8 %pi all 

(c226) mi : distrib (mi) ; 


f 2 (t ) + %i fl(t) f 2 (t ) %i fl(t) 

y ( + ) 

2 2 2 
(%i z - y) (- %i z - y) (- %i z - y) 



( d2 2 6 ) - 


8 %pi all 


%i fl(t) 


’si 


8 %pi all (%i z - y) 8 %pi all 


(c227) ki : first (mi) ; 


f 2 (t ) 


f 2 (t) + %i fl{t) 

y < + 

2 2 
(%i z - y) (- %i z - y) 


%i fl (t) 


(- %i z - y) 


( d2 2 7 ) 

(c228) ki : xthru (ki) ; 


8 %pi all 


(d228 ) - 


y ( ( f 2 (t ) - %i fl(t)) (%i z - y) + (f 2 (t) + %i fl (t) ) (- 

2 2 

8 %pi all (- %i z - y) (%i z - y) 


(c229) ki :expand(ki) ; 


( d2 2 9 ) 


2 f2 (t) y z 


4 2 2 4 

8 %pi all z +16 %pi all y z + 8 %pi all y 


4 fl(t) y z 


4 2 2 4 

8 %pi all z +16 %pi all y z + 8 %pi all y 


2 f 2 (t ) y 


4 2 2 4 

8 %pi all z +16 %pi all y z + 8 %pi all y 


(c230) ki : factor (ki) ; 


(d230) 


2 2 
y (f 2 (t) z +2 fl(t) y z - f2 (t) y ) 

2 2 2 

4 %pi all (z + y ) 

(c231) k j : rest (mi, 1) ; 

%i f'l(t) %i fl(t) 

( d2 3 1 ) 

8 %pi all (- %i z - y) 8 %pi all (%i z - y) 

(c232) k j : xthru (kj ) ; 

%i fl(t) (%i z - y) - %i fl(t) (- %i z - y) 

8 %pi all (- %i z - y) (%i z - y) 


(d232 ) 


fl(t) 


(- %i z - y) 


2 

%i z - y) ) 


(c233) kj :expand(kj) ; 



(d233 ) 

(c234) 

(d234) 

(c235) 

(d235) 

(c236) 

( d2 3 6 ) 

(c237 ) 
gtxyp: 

(d237 ) 

(c238 ) 

(d238 ) 

(c239) 

(d239) 


fl(t) z 


2 2 
4 %pi all z +4 %pi all y 

k j : factor (k j ) ; 

fl(t) z 


2 2 

4 %pi all (z + y ) 

mi : ki+k j ; 

2 2 

y (f2(t) z +2 fl(t) y z - f2(t) y ) fl(t) z 


2 2 2 2 2 
4 %pi all (z + y ) 4 %pi all (z + y ) 

mi : subst ( (t-x) ,z,mi) ; 

2 2 
y (- f2(t) y + 2 fl(t) (t-x) y + f2(t) (t-x) ) 


2 2 2 
4 %pi all (y + (t - x) ) 

fl(t) (t - x) 


2 2 
4 %pi all (y + (t - x) ) 

/* Get the answer of gsyyp expressed by the auxiliary function */ 
integrate (mi, t , -a, a) ; 
a 

/ 2 2 
[ y (- f 2 (t ) y + 2 fl(t) (t-x) y + f2 (t) (t-x) ) 

I ( 

] 2 2 2 
/ 4 %pi all (y + (t - x) ) 


fl(t) (t - x) 

) dt 

2 2 
4 %pi all (y + (t - x) ) 

ki:ki*4*%pi*all; 

2 2 
y (f 2 (t) z + 2 fl(t) y z - f2(t) y ) 


2 2 2 
(z + y ) 

kj :k j*4*%pi*all; 

fl(t) z 


2 2 
z + y 


(c240) mi : ki+k j; 



2 2 
y (f 2 (t ) z + 2 fl(t) y z - £2 (t) y ) fl(t) z 

( d2 4 0 ) 

2 2 2 2 2 
(z + y ) z + y 

(c241) gsxypf : subst ( (t-x) , z,mi) ; 

2 2 
y (- f 2 ( t ) y +2 fl(t) (t-x) y + f2 (t) (t-x) ) fl(t) (t - x) 

(d241) 

2 2 2 2 2 
(y + (t - x) ) y + (t - x) 

(c242 ) y : 0 ; 

(d242 ) 0 

(c243) ev(gtxyp); 

a 

/ 

( fl(t) 

I dt 

] t-x 

/ 

- a 

( d2 4 3 ) 

4 %pi all 

(c244 ) kill (y) ; 

(d244 ) done 

fl(t) : =1 ; 

(d2 45 ) f 1 (t) := 1 

(c246) f 2 (t ) : =0 ; 

(d246) f 2 (t ) := 0 

(c247) ''gsxxpf; 

2 2 

2 y y ( (t - x) - y ) 

( d2 4 7 ) + 

2 2 2 2 2 
y + (t - x) (y + (t - x) ) 

(c248) subst (yl , y, %) ; 

2 2 

2 yl yl ( (t - x) - yl ) 

( d2 4 8 ) + 

2 2 2 2 2 
yl + (t - x) (yl + (t - x) ) 

(c249) subst (xl, x, %) ; 

2 2 

2 yl yl ( (t - xl) - yl ) 

(d2 49) + 



2 


(c250 ) 
(d250) 
(c251) 
(d251) 
(c252) 
(d252 ) 
(c253) 


2 

yl + (t - xl) 
gsxxpf Itl : subst (tl, t , %) ; 


2 2 2 
(yl + (t - xl) ) 


2 2 

2 yl yl ((tl - xl) - yl ) 


2 2 2 2 2 
yl + (tl - xl) (yl + (tl - xl) ) 

' ' gsxxpf; 


2 2 

2 y y ( (t - x) - y ) 


2 2 2 2 2 
y + (t - x) (y + (t - x) ) 

subst (y2, y, %) ; 

2 2 

2 y2 y2 (<t - x) - y2 ) 


2 2 2 2 2 
y2 + (t - x) (y2 + (t - x) ) 

subst (x2, x, %) ; 


2 2 

2 y2 y2 ( (t - x2 ) - y2 ) 

(d253 ) + 

2 2 2 2 2 
y2 + (t - x2) (y2 + (t - x2) ) 

(c254) gsxxpflt2 : subst (t2,t, %) ; 


(d254 ) 

(c255) kill (f 1) ; 
(d255 ) 

(c256) kill (f 2) ; 
(d256) 

(c257) f 1 (t ) : =0 ; 
(d257) 

(c258) f2(t):=l; 
(d258) 

(c259) ' ' gsxxpf; 


2 2 

2 y2 y2 ( (t2 - x2) - y2 ) 


2 2 2 2 2 
y2 + (t2 - x2) (y2 + (t2 - x2) ) 


done 


done 


fl(t) := 0 
f 2 ( t ) := 1 


2 

t - x 2 (t - x) y 

(d259) 

2 2 2 2 2 
y + (t - x) (y + (t - x) ) 



(c260) subst (yl, y, %) ; 


2 


(d2 60 ) 

t - X 

2 (t 

- x) yl 

2 2 

2 

2 2 



yl + (t - x) 

(yl + 

(t - x) ) 

( c2 61 ) 

subst (xl, x, %) ; 


2 

(d26l) 

t - xl 

2 (t 

- xl) yl 

2 2 

2 

2 2 



yl + (t - xl) 

(yl + 

(t - xl) ) 

(c262 ) 

gsxxpf 2t 1 : subst (tl,t,%); 


2 

(d2 62 ) 

1 1 - xl 

2 (tl 

- xl) yl 

2 2 

2 

2 2 



yl + (tl - xl) 

(yl + 

(tl - xl) ) 

(c263) 

' ' gsxxpf; 


2 

(d263) 

t - X 

2 (t 

- x) y 

2 2 

2 

2 2 



y + (t - x) 

(y + (t - x) ) 

(c264) 

subst (y2, y, %) ; 


2 

(d2 64 ) 

t - X 

2 (t 

- x) y2 

2 2 

2 

2 2 



y2 + (t - x) 

(y2 + 

(t - x) ) 

(c265) 

subst (x2,x, %) ; 


2 

(d265) 

t — x2 

2 (t 

- x2 ) y2 

2 2 

2 

2 2 



y2 + (t - x2 ) 

(y2 + 

(t - x2) ) 

(c266) 

gsxxpf2t2 : subst (t2,t, %) ; 


2 

(d266) 

1 2 - x2 

2 (t2 

- x2) y2 

2 2 

2 

2 2 



y2 + (t2 - x2) 

(y2 + 

<t2 - x2 ) ) 

(c2 67 ) 

kill (fl) ; 



(d2 67 ) 


done 


(c268) 

kill (f 2) ; 




(d268) done 



fl(t) :=1; 


(d2 69) 
(c270 ) 
(d270) 
(c271) 

(d271) 

(c272 ) 

( 6212 ) 

(c273) 

( d2 7 3 ) 

(c274 ) 

(d274) 

(c275) 

(d275) 

(c276) 

(d276) 

(c277 ) 

( 6211 ) 


fl(t) := 1 

f 2 (t) : =0 ; 

f 2 (t) : = 0 

' ' gsyypf; 

2 2 

y ( (t - x) - y ) 


2 2 2 
(y + (t - x) ) 

subst (yl, y, %) ; 

2 2 
yl <(t - x) - yl ) 


2 2 2 
(yl + (t - x) ) 

subst (xl, x, %) ; 

2 2 
yl (<t - xl) - yl ) 


2 2 2 
(yl + <t - xl) ) 

gsyypfltl : subst (tl, t, %) ; 

2 2 

yl ((tl - xl) - yl ) 


2 2 2 
(yl + (tl - xl) ) 

' ' gsyypf; 

2 2 
y ( (t - x) - y ) 


2 2 2 
(y + (t - x) ) 

subst (y2,y, %) ; 

2 2 

y2 (<t - x) - y2 ) 


2 2 2 
(y2 + (t - x) ) 

subst (x2 , x, % ) ; 

2 2 

y2 ( (t - x2 ) - y2 ) 


2 2 2 
(y2 + (t - x2 ) ) 


(c278) gsyypf lt2 : subst (t2,t, %) ; 



(d278 ) 


2 2 

y2 ( (t2 - x2 ) - y2 ) 

2 2 2 
(y2 + (t2 - x2) ) 


(c279) kill (f 1) ; 
( d2 7 9 ) 

(c280 ) kill (f 2) ; 
( d2 8 0 ) 

(c281) fl(t):-0; 
( d2 81) 

(c282 ) f 2 (t ) : =1 ; 
(d282 ) 

(c283) ''gsyypf; 


done 


done 


fl(t) := 0 


f 2 (t ) := 1 


( d2 8 3 ) 


(c284) subst (yl,y,%) ; 


t - x 2 (t - x) y 

+ 

2 2 2 2 2 
y + (t - x) (y + (t - x) ) 


t - x 


( d2 8 4 ) 


2 (t - x) yl 


2 2 2 2 2 
yl + (t - x) (yl + (t - x) ) 


(c285) subst (xl, x, %) ; 


t - xl 


(d285) 


2 (t - xl) yl 


2 2 2 2 2 
yl + (t - xl) (yl + (t - xl) ) 


(c286) gsyypf 2tl : subst (tl, t, %} ; 


tl - xl 


( d2 8 6) 


2 (tl - xl) yl 


2 2 2 2 2 
yl + (tl - xl) (yl + (tl - xl) ) 


(c287 ) ' ' gsyypf; 


(d287 ) 


t - x 2 (t - x) y 

+ 

2 2 2 2 2 
y + (t - x) (y + (t - x) ) 


(c288) subst (y2 , y, %) ; 


2 



(d288 ) 

t - X 


2 

(t 

- x) y2 

2 2 


2 


2 2 



y2 + (t - x) 


(y2 

+ 

(t - x) ) 

(c289) 

subst (x2 , x, %) ; 




2 

(d289 ) 

t - x2 


2 

(t 

- x2 ) y2 

2 2 


2 


2 2 



y2 + (t - x2) 


(y2 

+ 

(t - x2 ) ) 

(c290) 

gsyypf2t2 : subst (t2,t, %) ; 




2 

(d290 ) 

1 2 - x2 


2 

(t2 

- x2) y2 

2 2 


2 


2 2 



y2 + (t2 - x2) 


(y2 

+ 

(t2 - x2 ) ) 


(c291) kill (f 1) ; 

(d2 91 ) 

(c292) kill (f 2) ; 

( d2 92) 

(c293) 
fl(t) :=1; 

(d2 93) 

(c294 ) f2(t):=0; 

( d2 9 4 ) 

(c295) ''gsxypf; 

(d2 95 ) 

(c296) subst (yl, y, %) ; 
( d2 9 6 ) 

(c297) subst (xl, x, %) ; 
(d297 ) 


done 


done 



fl(t) 

:= 1 


f 2 (t ) 

:= 0 


2 


2 (t - x) 

y 

t - X 

2 

2 2 

2 2 

(y + (t - 

x) ) 

y + (t - x) 


2 


2 (t - x) 

yi 

t - X 

2 

2 2 

2 2 

(yl + (t - 

x) ) 

yl + (t - x) 


2 


2 (t - xl) 

yi 

t - xl 


2 2 2 2 2 
(yl + (t - xl) ) yl + (t - xl) 


(c298) gsxypf ltl : subst (tl , t, %) 


2 



(d2 98 ) 
(c299) 
(d299) 
(c300 ) 
( d3 0 0 ) 
(c301) 
(d301) 
(c302 ) 
(d302 ) 


2 (tl - xl) yl tl - xl 


2 2 2 2 2 
(yl + (tl - xl) ) yl + (tl - xl) 

' ' gsxypf; 

2 

2 (t - x) y t - x 


2 2 2 2 2 
(y + (t - x) ) y + (t - x) 

subst (y2, y, %) ; 

2 

2 (t - x) y2 t - x 


2 2 2 2 2 
(y2 + (t - x) ) y2 + (t - x) 

subst (x2, x, %) ; 

2 

2 (t - x2) y2 t - x2 


2 2 2 2 2 
(y2 + (t - x2) ) y2 + (t - x2) 

gsxypf It 2 : subst (t2 , t , %) ; 

2 

2 ( 1 2 - x2 ) y2 t2 - x2 


2 2 2 2 2 
(y2 + (t2 - x2) ) y2 + (t2 - x2) 


(c303) kill (f 1) ; 
(d303) 

(c304 ) kill (f 2) ; 

(d304 ) 

(c305) fl (t) :=0; 

(d305 ) 

(c306) f 2 (t ) : =1 ; 

( d3 0 6 ) 

(c307) '' gsxypf; 

(d307 ) 

(c308) subst (yl, y, %) ; 


done 


done 


fl(t) := 0 


f 2 (t) := 1 


2 2 
y (<t - x) - y ) 


2 2 2 
(y + (t - x) ) 


2 2 
yl ( (t - x) - yl ) 


( d3 0 8 ) 



(c309) subst (xl, x, %) ; 


2 2 2 
(yl + (t - x) ) 


(d309) 


2 2 

yl ((t - xl) - yl ) 

2 2 2 
(yl + (t - xl) ) 


(c310) gsxypf2tl : subst (tl, t, %) ; 


(d310 ) 


(c311) ''gsxypf; 


(d311) 


(c312) subst (y2,y,%); 


(d312) 


(c313) subst (x2, x, %) ; 


(d313) 


2 2 

yl ((tl - xl) - yl ) 

2 2 2 
(yl + (tl - xl) ) 


2 2 
y ( (t - x) - y ) 

2 2 2 
(y + (t - x) ) 


2 2 

y2 ((t - x) - y2 ) 

2 2 2 
(y2 + (t - x) ) 


2 2 
y2 ( (t - x2) - y2 ) 

2 2 2 
(y2 + (t - x2 ) ) 


(c3l4) gsxypf2t2 : subst (t2,t,%); 


(d314 ; 


2 2 

y2 ( (t2 - x2) - y2 ) 

2 2 2 
(y2 + (t2 - x2) ) 


(c315 ) kill (f 1) ; 

(d3l5) done 

(c316) kill (f 2) ; 

(d316) done 

(c317) /***** calculate kernel of two cracks *************************************** 
kill ( wwm) ; 

(d317) done 

(c318) xl:al*gc; 



( d3 18) 

(c319) tl:al*gt; 

( d3 1 9 ) 

(c320) t2:a2*gt; 

(d320 ) 

(c321) gsx2 : ' ' gsxxpf lt2 ; 


al gc 


al gt 


a2 gt 


(d321) 


2 y2 


2 2 
y2 + (a2 gt - x2) 


2 2 

y2 ( (a2 gt - x2) - y2 ) 

2 2 2 
(y2 + (a2 gt - x2) ) 


(c322) gsy2 : ' ' gsyypf lt2 ; 


(d322 ) 


2 2 

y2 ( (a2 gt - x2) - y2 ) 

2 2 2 
(y2 + (a2 gt - x2) ) 


(c323) gtxy2 : ' ' gsxypf lt2; 


(d323) 


2 (a2 gt - x2) y2 


a2 gt - x2 


2 2 2 2 2 
(y2 + (a2 gt - x2) ) y2 + (a2 gt - x2) 


(c324) x2 :p3+xl*cos (-gh) -yl*sin (-gh) ; 


(d324) 


sin(gh) yl + p3 + al gc cos (gh) 


(c325) y2 :p4+xl*sin (-gh) +yl*cos (-gh) ; 


( d3 2 5 ) 


cos (gh) yl + p4 - al gc sin(gh) 


(c326) /* 1 * k(ll,21) 2 to 1 term f21 gtxy *** c21f21xy *************************** 
gtxy21b : - ( ' gsx2- f gsy2 ) /2*sin ( 2* (-gh) ) +' gtxy 2* cos (2* (-gh) ) ; 


(d326) 


cos (2 gh) gtxy2 + 


sin(2 gh) (gsx2 - gsy2) 


( c327 ) gtxy 21 : - ( ' ' gsx2- f 'gsy2) /2*sin(2* (-gh) ) +' ' gtxy2*cos (2* ( -gh) ) ; 


(d327) cos(2 gh) (2 (cos(gh) yl + p4 - al gc sin(gh)) 


(- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 


/expt((- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 


+ (cos(gh) yl + p4 - al gc sin(gh)) , 2) 


(- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 



2 


/((- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

+ (cos(gh) yl + p4 - al gc sin(gh)) )) 

+ sin(2 gh) (2 (cos(gh) yl + p4 - al gc sin(gh)) 

2 

/((- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

+ (cos (gh) yl + p4 - al gc sin(gh)) ) + 2 (cos (gh) yl + p4 

2 

((- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 


- (cos (gh) yl + p4 - al gc sin(gh)) ) 

2 

/expt((- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

+ (cos (gh) yl + p4 - al gc sin (gh) ) , 2))/2 
(c328) yl : 0 ; 

( d3 28) 0 

(c329) ev(gtxy21); 

2 (p4 - al gc sin(gh)) 


(d329) sin (2 gh) (- 


(p4 - al gc sin(gh)) + (- p3 + a2 gt - al 

2 

+ 2 (p4 - al gc sin(gh)) ( (- p3 + a2 gt - al gc cos(gh)) 


- (p4 - al gc sin(gh)) ) 

2 2 2 
/ ( (p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos (gh) ) ) )/ 


+ cos (2 gh) (- 


2 (- p3 + a2 gt - al gc cos (gh) ) (p4 - al gc 

2 

( (p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc 
- p3 + a2 gt - al gc cos(gh) 


•) 


2 2 
(p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos(gh)) 

(c330) subst (p3a2al, -p3+a2*gt-al*gc*cos (gh) , %) ; 

2 (p4 - al gc sin(gh)) 


(d330) sin (2 gh) (- 


2 2 
(p4 - al gc sin(gh)) + p3a2al 


2 (p4 - al gc sin(gh)) (p3a2al 


- al gc sin (gh) ) 


2 

gc cos (gh) ) 


2 

2 

sin (gh) ) 


2 2 

cos (gh) ) ) 


2 


2 

(p4 - al gc sin(gh)) ) 



2 2 2 

( (p4 - al gc sin(gh)) + p3a2al ) 


+ cos (2 gh) (- 


2 p3a2al (p4 - al gc sin(gh)) 

2 2 2 

( (p4 - al gc sin(gh)) + p3a2al ) 
p3a2al 


2 2 

(p4 - al gc sin(gh)) + p3a2al 
(c331) subst (p4al,p4-al*gc*sin (gh) , %) ; 


sin(2 gh) ( 


(d331 ) 


2 2 

2 p4al 2 p4al (p3a2al - p4al ) 

+ ) 

2 2 2 2 2 
p4al + p3a2al (p4al + p3a2al ) 


2 p3a2al p4al p3a2al 

+ cos (2 gh) ( ) 

2 2 2 2 2 
(p4al + p3a2al ) p4al + p3a2al 


(c332) subst (ww,p4al A 2+p3a2al A 2, %) ; 

2 2 

2 p4al 2 p4al (p3a2al - p4al ) 

sin(2 gh) ( + ) 

ww 2 


(d332 ) 


2 p3a2al p4al p3a2al 

+ cos (2 gh) ( ) 

2 ww 


(c3 33) fn : subst (wwm, p3a2al A 2-p4al A 2, %) ; 

2 p4al wwm 2 p4al 

sin(2 gh) ( + ) 

2 ww 


(d333) 


2 p3a2al p4al p3a2al 

+ cos (2 gh) ( ) 

2 ww 


(c334) first (fn) ; 


2 p4al wwm 2 p4al 

sin (2 gh) ( + ) 

2 ww 


(d334) 


2 



(c335) n2 :multthru {%) ; 


sin (2 gh) p4al wwm sin (2 gh) p4al 


(d335) + 

2 ww 

ww 


(c336) rest(fn,l); 

2 


2 p3a2al p4al p3a2al 

(d336) cos (2 gh) ( ) 

2 ww 


ww 


(c337) n3 :multthru (%) ; 


(d337 ) 


2 

2 cos (2 gh) p3a2al p4al 


2 

ww 


cos (2 gh) p3a2al 


ww 


(c338 ) 
(d338 ) 


n : n2+n3; 

sin (2 gh) p4al wwm sin (2 gh) p4al 

+ 

2 ww 

ww 


cos (2 gh) p3a2al 


ww 


2 

2 cos (2 gh) p3a2al p4al 

+ 

2 

ww 

(c339) fn : combine (n) ; 

2 

sin (2 gh) p4al wwm + 2 cos (2 gh) p3a2al p4al 

(d339) 

2 

ww 


sin (2 gh) p4al - cos (2 gh) p3a2al 
+ 


(c340) /* 1 * k (11, 21) 2 to 1 term f21 gtxy *** 
kll21: (a2/ (4*%pi*all) )*fn; 


ww 


2 

sin(2 gh) p4al wwm + 2 cos (2 gh) p3a2al p4al 

(d340 ) a2 ( 

2 

ww 

sin{2 gh) p4al - cos (2 gh) p3a2al 

+ ) / (4 %pi all) 

ww 


(c341) /* where 

p3a2al=-p3+a2*gt-al*gc*cos (gh) 


p4al=p4-al*gc*sin (gh) 



ww=p4al A 2+p3a2al^2 

wwm=p3a2al A 2-p4al / '2 

*/ 

fortran (p3a2al : -p3+a2*gt ( j ) -al*gc (i) *dcos (gh) ) $ 
p3a2al = -p3+a2*gt ( j) -al*dcos (gh) *gc (i) 

(c342) fortran (p4al :p4-al*gc (i) *dsin (gh) ) $ 
p4al = p4-al*dsin (gh) *gc (i) 

(c343) fortran (ww : ' p4al^2+' p3a2al A 2 ) $ 
ww = p4al**2+p3a2al**2 

(c34 4 ) fortran (wwm: 'p3a2al A 2-' p4al A 2) $ 
wwm = p3a2al**2-p4al**2 

(c345) subst (dcos, cos, first (first (fn) ) ) $ 

(c346) subst (dsin, sin, %) $ 

(c347) fortran (part 1 :%) $ 

parti = dsin (2*gh) *p4al*wwm+2*dcos (2*gh) *p3a2al*p4al**2 

(c348 ) fortran (parti : ' parti/ ( ' ww A 2) ) $ 
parti = partl/ww**2 

(c349) subst (dcos, cos, last (fn) ) $ 

(c350) subst (dsin, sin, %) $ 

(c351) fortran (part2 :%) $ 

part2 = (dsin (2*gh) *p4al-dcos (2*gh) *p3a2al) /ww 

( c352 ) kill (xl) $ 

( c353 ) kill (x2) $ 

(c354 ) kill (yl) $ 

(c355) kill (y2 ) $ 

(c356) kill (tl) $ 

(c357) kill (t2) $ 

(c358) kill (p3a2al) $ 

(c359) kill(p4al)$ 

(c360) kill (ww) $ 

(c361) kill (wwm) $ 

(c362) xl:al*gc; 

(d362) al gc 

(c363) tl:al*gt; 

( d3 63) al gt 

(c364) t2 : a2*gt ; 


(d364) 


a2 gt 



(c365) gsx2 : ' ' gsxxpf 2t2 ; 


(d365) 


(c366) 


a2 gt - x2 

2 2 
y2 + (a2 gt - x2) 

gsy2 : ' ' gsyypf 2t2 ; 


2 

2 (a2 gt - x2) y2 

2 2 2 
(y2 + (a2 gt - x2) ) 


( d3 6 6 ) 


a2 gt - x2 


2 (a2 gt - x2) y2 


2 2 2 2 2 
y2 + {a2 gt - x2) (y2 + (a2 gt - x2) ) 


(c367) gtxy2 : ' ' gsxypf 2t2 ; 


2 2 


(d3 67 ) 


y2 ( (a2 gt - x2) - y2 ) 

2 2 2 


(y2 + (a2 gt - x2) ) 

(c368) x2 :p3+xl*cos (-gh) -yl*sin (-gh) ; 

(d368) sin(gh) yl + p3 + al gc cos (gh) 

(c369) y2 :p4+xl*sin (-gh) +yl*cos (-gh) ; 

(d369) cos (gh) yl + p4 - al gc sin(gh) 

(c370 ) /* 2 * k (11, 22) 2 to 1 term £22 gtxy *** c21f22xy *************************** 

gtxy21b : - ( ' gsx2-' gsy2 ) /2*sin (2* (-gh) ) +' gtxy2*cos (2* (-gh) ) ; 

sin(2 gh) (gsx2 - gsy2) 

(d370) cos (2 gh) gtxy2 + 


2 

(c371) gtxy 21 (' ' gsx2-' 'gsy2)/2*sin(2*(-gh))+' 'gtxy2*cos (2* (-gh) ) ; 
(d371) cos (2 gh) (cos (gh) yl + p4 - al gc sin(gh)) 


2 

( (- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

- (cos (gh) yl + p4 - al gc sin(gh)) ) 

2 

/expt ( (- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

+ (cos(gh) yl + p4 - al gc sin(gh)) , 2) 

2 

- 2 sin (2 gh) (cos (gh) yl + p4 - al gc sin (gh) ) 

(- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 

/ expt ( (- sin (gh) yl - p3 + a2 gt - al gc cos(gh)) 


2 



+ (cos (gh) yl + p4 - al gc sin(gh)) , 2) 

(c372 ) yl : 0 ; 

(d372 ) 0 

(c373) ev(gtxy21); 

2 

(d373) cos{2 gh) (p4 - al gc sin(gh)) ( (- p3 + a2 gt - al gc cos(gh)) 

2 


- (p4 - al gc sin(gh)) ) 

2 2 2 
/ ( (p4 - al gc sin (gh) ) + (- p3 + a2 gt - al gc cos (gh) ) ) 

2 

2 sin (2 gh) (- p3 + a2 gt - al gc cos(gh)) (p4 - al gc sin(gh)) 

2 2 2 


(c374) 


( d3 7 4 ) 


( (p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos(gh)) ) 

subst (p3a2al, -p3+a2*gt-al*gc*cos (gh) , %) ; 

2 2 
cos (2 gh) (p4 - al gc sin(gh)) (p3a2al - (p4 - al gc sin(gh)) ) 


2 2 2 

( (p4 - al gc sin(gh)) + p3a2al ) 

2 

2 sin (2 gh) p3a2al (p4 - al gc sin(gh)) 


2 2 2 

( (p4 - al gc sin(gh)) + p3a2al ) 


(c375) subst (p4al,p4-al*gc*sin (gh) , %) ; 


(d375) 


2 2 2 
cos (2 gh) p4al (p3a2al - p4al ) 2 sin (2 gh) p3a2al p4al 


2 2 2 
(p4al + p3a2al ) 


2 2 2 
(p4al + p3a2al ) 


(c376) subst (ww,p4al A 2+p3a2al A 2 , %) ; 


(d376) 


2 2 2 
cos (2 gh) p4al (p3a2al - p4al ) 2 sin (2 gh) p3a2al p4al 


ww 


ww 


(c377) f n : subst (wwm, p3a2al A 2-p4al A 2 , %) ; 


(d377 ; 


cos (2 gh) p4al wwm 2 sin (2 gh) p3a2al p4al 


ww 


ww 


(c378) fn : combine (fn) ; 


cos (2 gh) p4al wwm - 2 sin (2 gh) p3a2al p4al 


( d3 7 8 ) 



2 


ww 

(c379) /* 2 * k(ll,22) 2 to 1 term f22 gtxy *** c21f22xy *************************** 
kl!22 : (a2/ (4*%pi*all) ) *fn; 


( d3 7 9 ) 


2 

a2 (cos (2 gh) p4al wwm - 2 sin(2 gh) p3a2al p4al ) 


2 

4 %pi all ww 


(c380) /* where 

p3a2al=-p3+a2*gt-al*gc*cos (gh) 

p4al=p4-al*gc*sin (gh) 

ww=p4al A 2+p3a2al A 2 

wwm=p3a2al A 2-p4al A 2 

*/ 

length (fn) ; 

(d380 ) 2 

(c381) fortran (p3a2al : -p3+a2*gt ( j) -al*gc (i) *dcos (gh) ) $ 
p3a2al = -p3+a2*gt ( j) -al*dcos (gh) *gc (i) 

(c382) fortran (p4al :p4-al*gc (i) *dsin (gh) ) $ 
p4al = p4-al*dsin (gh) *gc (i) 

(c383) fortran (ww : ' p4al A 2+' p3a2al A 2 ) $ 
ww = p4al**2+p3a2al**2 

(c384) fortran (wwm: / p3a2al A 2- , p4al A 2) $ 
wwm = p3a2al**2-p4al**2 

(c385) subst (dcos, cos , f irst (f n) ) $ 

(c386) subst (dsin, sin, %) $ 

(c387) fortran (part 1 :%) $ 

parti = dcos (2*gh) *p4al*wwm-2*dsin (2*gh) *p3a2al*p4al**2 

(c388) fortran (parti : 'parti/ ( ' ww A 2) ) $ 
parti = partl/ww**2 

(c389) kill (xl) $ 

(c390) kill (x2 ) $ 

(c391) kill (yl) $ 

(c392 ) kill (y2) $ 

(c393) kill (tl) $ 

(c394 ) kill (t2 ) $ 

(c395) kill (p2a2al) $ 

(c396) kill (p4al) $ 

(c397) kill (ww) $ 





(c398) kill (wwm) $ 

(c399) xl:al*gc; 

( d3 9 9 ) al gc 

(c400) tl:al*gt; 

(d400 ) al gt 

(c401) t2:a2*gt; 

(d401) a2 gt 

(c402) gsx2 : ' ' gsxxpf lt2 ; 

2 2 

2 y2 y2 ( (a2 gt - x2) - y2 ) 

(d402 ) + 

2 2 2 2 2 
y2 + (a2 gt - x2) (y2 + (a2 gt - x2) ) 

(c403) gsy2 : ' 'gsyypflt2; 


(d403) 


2 2 

y2 ( (a2 gt - x2) - y2 ) 


2 2 2 
(y2 + (a2 gt - x2) ) 


(c404 ) gtxy2 : ' ' gsxypf lt2; 


(d404 ) 


2 (a2 gt - x2) y2 


a2 gt - x2 


2 2 2 2 2 
(y2 + (a2 gt - x2) ) y2 + (a2 gt - x2) 


(c405) x2 :p3+xl*cos (-gh) -yl*sin (-gh) ; 

(d405) sin (gh) yl + p3 + al gc cos (gh) 

(c406) y2 :p4+xl*sin (-gh) +yl*cos (-gh) ; 

(d 406) cos (gh) yl + p4 - al gc sin (gh) 

(c407) /* 3 * k(12,21) 2 to 1 term f21 gsyy *** c21f21y ************************** 
gsy21b : (' gsx2+' gsy2 ) / 2- ( ' gsx2-' gsy2) /2*cos (2* (-gh) ) -' gtxy2*sin (2* (-gh) ) ; 

gsy2 + gsx2 cos (2 gh) (gsx2 - gsy2) 

2 2 
(c4 08 ) gsy21 : ( ' ' gsx2+' ' gsy2)/2-(''gsx2-''gsy2)/2*cos(2* (-gh) ) -' ' gtxy2*sin (2 * (-gh) ) 

2 

(d408) sin(2 gh) (2 (cos(gh) yl + p4 - al gc sin(gh)) 

(- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 

/expt ( (- sin (gh) yl - p3 + a2 gt - al gc cos (gh) ) 


(d407 ) 


sin (2 gh) gtxy2 + 


2 



+ (cos(gh) yl + p4 - al gc sin(gh)) , 2) 


- (- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 

/((- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 

+ (cos (gh) yl + p4 - al gc sin(gh)) )) 

- cos (2 gh) (2 (cos (gh) yl + p4 - al gc sin(gh)) 

2 

/ ( (- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

+ (cos(gh) yl + p4 - al gc sin(gh)) ) + 2 (cos (gh) yl + p4 - al gc sin(gh)) 

2 

((- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 

- (cos (gh) yl + p4 - al gc sin(gh)) ) 

2 

/expt((- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 

+ (cos (gh) yl + p4 - al gc sin(gh)) , 2 ) ) / 2 
+ (cos (gh) yl + p4 - al gc sin(gh)) 

2 

/ ( (- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

+ (cos (gh) yl + p4 - al gc sin(gh)) ) 

(c409) yl:0; 

(d409) 0 

(c410) ev(gsy21); 

2 (p4 - al gc sin(gh)) 

(d410) - cos (2 gh) ( 

2 2 
(p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos (gh) ) 

2 

+ 2 (p4 - al gc sin(gh)) ( (- p3 + a2 gt - al gc cos(gh)) 

2 

- (p4 - al gc sin(gh)) ) 

2 2 2 
/ ( (p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos (gh) ) ) )/2 


+ sin(2 gh) (■ 


z 

2 (- p3 + a2 gt - al gc cos(gh)) (p4 - al gc sin(gh)) 

2 2 2 
( (p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos(gh)) ) 

- p3 + a2 gt - al gc cos(gh) 
) 



2 


2 

(p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos(gh)} 

p4 - al gc sin(gh) 

+ 

2 2 
(p4 - al gc sin(gh)) + <- p3 + a2 gt - al gc cos(gh)) 

(c411) subst (p3a2al, -p3+a2*gt-al*gc*cos (gh) , %) ; 

(d411) - cos (2 gh) 


(- 


2 (p4 - al gc sin(gh)) 


2 2 
(p4 - al gc sin(gh)) + (- p3 + a2 gt(j) - al dcos(gh) gc(i)) 


+ 2 (p4 - al gc sin(gh)) ( (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 

2 

- (p4 - al gc sin(gh)) ) 

2 2 2 
/ ( (p4 - al gc sin(gh)) + (- p3 + a2 gt{j) - al dcos(gh) gc(i)) ) )/2 


+ sin(2 gh) (- 


2 (- p3 + a2 gt(j) - al dcos(gh) gc(i)) (p4 - al gc sin(gh)) 

2 2 2 
( (p4 - al gc sin(gh)) + {- p3 + a2 gt(j) - al dcos (gh) gc(i)) ) 


- p3 + a2 gt(j) - al dcos (gh) gc(i) 


2 2 
(p4 - al gc sin(gh)) + (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 

p4 - al gc sin(gh) 


-) 


+ 


2 2 
(p4 - al gc sin(gh)) + (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 

(c412) subst (p4al,p4-al*gc*sin (gh) , %) ; 


2 p4al 

(d412 ) - cos (2 gh) ( 

2 2 
p4al + (- p3 + a2 gt(j) - al dcos(gh) gc(i)) 

2 2 

2 p4al ( (- p3 + a2 gt(j) - al dcos (gh) gc(i)) - p4al ) 

+ ) / 2 

2 2 2 


(p4al + (- p3 + a2 gt(j) - al dcos (gh) gc(i)) ) 

2 


+ sin(2 gh) 


2 (- p3 + a2 gt(j) - al dcos(gh) gc(i)) p4al 


2 


2 2 


(p4al + (- p3 + a2 gt(j) - al dcos(gh) gc(i)) ) 


- p3 + a2 gt(j) - al dcos (gh) gc(i) 



p4al + (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 



p4al 



p4al + (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 


(c413) subst (ww, p4al A 2+p3a2al A 2 , %) ; 
(d4 13 ) 


cos (2 gh) 


2 2 

2 p4al 2 p4al ((- p3 + a2 gt(j) - al dcos (gh) gc(i)) - p4al ) 

( + 

ww 2 


ww 


2 


+ sin(2 gh) 


2 

2 (- p3 + a2 gt(j) - al dcos(gh) gc(i)) p4al 


ww 


- p3 + a2 gt(j) - al dcos(gh) gc(i) p4al 

) + 

ww ww 

( c414 ) f n : subst (wwm, p3a2al A 2-p4al A 2 , %) ; 

2 p4al wwm 2 p4al 

cos (2 gh) ( + ) 

2 ww 


(d414) 


ww 


2 


+ sin (2 gh) 


2 

2 (- p3 + a2 gt(j) - al dcos(gh) gc(i)) p4al 

( 

2 


ww 

- p3 + a2 gt(j) - al dcos (gh) gc(i) p4al 

) + 

ww ww 


(c415) nl:last(fn); 
(d415) 

(c416) first (fn) ; 


p4al 

ww 


( d4 1 6 ) 

(c417) n2 :multthru (%) ; 


2 p4al wwm 2 p4al 

cos (2 gh) ( + ) 

2 ww 

ww 


2 


cos (2 gh) p4al wwm cos (2 gh) p4al 


(d4 17 ) 


2 


ww 



ww 


(c418) rest {fn, 1) ; 


(d418) sin(2 gh) 


2 

2 (- p3 + a2 gt(j> - al dcos (gh) gc(i)) p4al 

( 

2 

ww 

- p3 + a2 gt ( j) - al dcos (gh) gc(i) p4al 

) + 

ww ww 


(c419) first (%) ; 


(d419) sin(2 gh) 


2 

2 (- p3 + a2 gt ( j) - al dcos(gh) gc(i)) p4al 

( 

2 

ww 


- p3 + a2 gt(j) - al dcos (gh) gc(i) 

} 

ww 


(c420) n3 : mult thru (%) ; 


2 

2 sin(2 gh) (- p3 + a2 gt(j) - al dcos(gh) gc(i)) p4al 

(d420 ) 

2 

ww 

sin (2 gh) (- p3 + a2 gt(j) - al dcos(gh) gc(i)) 


ww 


(c421) n:nl+n2+n3; 


cos (2 gh) p4al wwm cos (2 gh) p4al p4al 

(d4 21) + 

2 ww ww 

ww 


sin (2 gh) (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 


ww 


2 

2 sin (2 gh) (- p3 + a2 gt(j) - al dcos (gh) gc(i)) p4al 

+ 

2 

ww 


(c422) n : combine (n) / 

(d422 ) 

2 

2 sin (2 gh) (- p3 + a2 gt(j) - al dcos(gh) gc(i)) p4al - cos (2 gh) p4al wwm 


2 

ww 

- cos (2 gh) p4al + p4al - sin(2 gh) (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 



+ 


ww 

(c423) tl : first (% ) ; 

(d423) 

2 

2 sin(2 gh) (- p3 + a2 gt(j) - al dcos (gh) gc(i)) p4al - cos(2 gh) p4al wwm 


2 

ww 

(c424) last (n) ; 

(d424 ) 

- cos (2 gh) p4al + p4al - sin (2 gh) (- p3 + a2 gt(j) - al dcos (gh) gc(i)) 


ww 

(c425) t2 :part (%, 1) ; 

(d425) - cos (2 gh) p4al + p4al - sin(2 gh) 

(- p3 + a2 gt ( j) - al dcos(gh) gc(i)) 

(c426) first (%) +first (rest (%, 1) ) ; 

(d426) p4al - cos (2 gh) p4al 

( c427 ) factor (%) ; 

(d427 ) - (cos (2 gh) - 1) p4al 

(c428) trigexpand (% ) ; 

2 2 

(d428) - (- sin (gh) + cos (gh) - 1) p4al 

(c429) t2 : (trigsimp (%) +last ( 1 2 ) ) /ww; 

2 

2 sin (gh) p4al - sin(2 gh) (- p3 + a2 gt(j) - al dcos (gh) gc(i)} 

(d429) 

ww 

(c430) /* 3 * k(12,21) 2 to 1 term f21 gsyy *** c21f21y **************************** 
kl221 : (a2/ (4*%pi*all) ) * (tl+t2) ; 

(d430 ) a2 

2 

2 sin (2 gh) (- p3 + a2 gt(j) - al dcos(gh) gc(i)) p4al - cos (2 gh) p4al wwm 

( 

2 

ww 


2 

2 sin (gh) p4al - sin (2 gh) (- p3 + a2 gt ( j ) - al dcos (gh) gc(i)) 

+ ) 

ww 


/ (4 %pi all) 



(c431) /* where 

p3a2al=-p3+a2*gt-al*gc*cos (gh) 

p4al=p4-al*gc*sin (gh) 

ww=p4al A 2+p3a2al A 2 

wwm=p3a2al A 2-p4al A 2 
*/ 

fortran (p3a2al : -p3+a2*gt ( j ) -al*gc (i) *dcos (gh) ) $ 
p3a2al = -p3+a2*gt ( j) -al*dcos (gh) *gc (i) 

(c432) fortran (p4al :p4-al*gc (i) *dsin (gh) ) $ 
p4al = p4-al*dsin (gh) *gc (i) 

(c433 ) fortran (ww : ' p4al A 2+' p3a2al A 2 ) $ 
ww = p4al**2+p3a2al**2 

(c434) fortran (wwm: 'p3a2al A 2-'p4al A 2) $ 
wwm = p3a2al**2-p4al**2 

(c435) subst (dcos, cos, first (tl) ) $ 

(c436) subst (dsin, sin, %) $ 

(c437) fortran (parti :%) $ 

parti = 2*dsin (2*gh) * (-p3+a2*gt ( j ) -al*dcos (gh) *gc (i) ) *p4al**2-dcos 
1 (2*gh) *p4al*wwm 

(c438 ) fortran (parti : 'parti/ ( ' ww A 2) ) $ 
parti = partl/ww**2 

(c439) subst (dcos , cos , t2) $ 

(c440) subst (dsin, sin, %) $ 

(c441) fortran (part2 :%) $ 

part 2 = (2*dsin (gh) **2*p4al-dsin (2*gh) * (-p3+a2*gt ( j ) -al*dcos (gh) *g 
1 c (i) ) ) /ww 

(c442) kill (xl ) $ 

(c443) kill (x2) $ 

(c444) kill (yl) $ 

(c445) kill (y2) $ 

(c446) kill (t 1) $ 

(c447) kill (t2) $ 

(c448 ) kill (p3a2al) $ 

(c449) kill (p4al) $ 

(c450) kill (ww) $ 

(c451) kill (wwm) $ 

(c452) xl:al*gc; 

(d452) al gc 


(c453) tl:al*gt; 



(d4 53 ) 


al gt 


(c454) t2:a2*gt; 

(d454) a2 gt 

(c455) gsx2 : ' ' gsxxpf 2t2 ; 


2 

a2 gt - x2 2 (a2 gt - x2) y2 

(d 4 55) 

2 2 2 2 2 

y2 + (a2 gt - x2) (y2 + (a2 gt - x2) ) 

(c456) gsy2 : ' ' gsyypf 2t2 ; 

2 

a2 gt - x2 2 (a2 gt - x2) y2 

(d4 56 ) + 

2 2 2 2 2 

y2 + (a2 gt - x2) (y2 + (a2 gt - x2) ) 

(c457) gtxy2 : ' ' gsxypf 2t2; 


(d457 ) 

(c458 ) 
(d458) 
(c459) 
(d459) 
(c460) 
gsy2lb 

(d4 60 ) 


2 2 

y2 ( (a2 gt - x2) - y2 ) 

2 2 2 
(y2 + (a2 gt - x2) ) 

x2 :p3+xl*cos (-gh) -yl*sin (-gh) ; 

sin(gh) yl + p3 + al gc cos (gh) 

y2 :p4+xl*sin (-gh) +yl*cos (-gh) ; 

cos (gh) yl + p4 - al gc sin (gh) 

/* 4 * k(12,22) 2 to 1 term f22 gsyy *** c21f22y ************************** 

( ' gsx2 + ' gsy2 ) /2- ( ' gsx2-' gsy2 ) /2*cos (2* (-gh) ) -' gtxy2*sin (2* (-gh) ) ; 

gsy2 + gsx2 cos (2 gh) (gsx2 - gsy2) 

sin (2 gh) gtxy2 + 

2 2 


( c4 61 ) gsy21 : ( ' ' gsx2+' , gsy2)/2-( ,, gsx2- ,, gsy2)/2*cos(2* (-gh) ) -' ' gtxy2*sin (2 * ( -gh) ) 
(d46l) (- sin (gh) yl - p3 + a2 gt - al gc cos (gh) ) 

2 

/ ( (- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

+ (cos (gh) yl + p4 - al gc sin(gh)) ) + sin (2 gh) 

(cos (gh) yl + p4 - al gc sin(gh)) 


2 

((- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 

2 

- (cos (gh) yl + p4 - al gc sin(gh)) ) 

2 



/expt((- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 


+ (cos (gh) yl + p4 - al gc sin(gh)) , 2) 

2 

+ 2 cos (2 gh) (cos (gh) yl + p4 - al gc sin(gh)) 
(- sin(gh) yl - p3 + a2 gt - al gc cos (gh) ) 

/expt((- sin(gh) yl - p3 + a2 gt - al gc cos(gh)) 


+ (cos (gh) yl + p4 - al gc sin(gh)) , 2) 

(c462) yl:0; 

(d4 62 ) 0 

(c463) ev (gsy21) ; 

- p3 + a2 gt - al gc cos ( gh) 


(d463) 


2 2 
(p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos(gh)) 


+ sin (2 gh) (p4 - al gc sin(gh)) ( (- p3 + a2 gt - al gc cos(gh) 


- (p4 - al gc sin(gh)) ) 

2 2 2 
/ ( (p4 - al gc sin (gh) ) + (- p3 + a2 gt - al gc cos (gh) ) ) 


2 cos (2 gh) (- p3 + a2 gt - al gc cos(gh)) (p4 - al gc sin(gh)) 

+ 

2 2 2 
( (p4 - al gc sin(gh)) + (- p3 + a2 gt - al gc cos(gh)) ) 

(c464) subst (p3a2al, -p3+a2*gt-al*gc*cos (gh) , %) ; 

p3a2al 


(d464 ) 


2 2 

(p4 - al gc sin(gh)) + p3a2al 


2 2 
sin (2 gh) (p4 - al gc sin(gh)) (p3a2al - (p4 - al gc sin(gh)) ) 

2 2 2 

( (p4 - al gc sin(gh)) + p3a2al ) 


2 cos (2 gh) p3a2al (p4 - al gc sin(gh)) 


+ 


2 2 2 

( (p4 - al gc sin(gh)) + p3a2al ) 
(c4 65) subst (p4al,p4-al*gc*sin (gh) , %) ; 


(d4 65) 


p3a2al 


2 cos (2 gh) p3a2al p4al 


+ 



2 2 2 2 2 
p4al + p3a2al (p4al + p3a2al ) 

2 2 

sin (2 gh) p4al (p3a2al - p4al ) 


2 2 2 
(p4al + p3a2al ) 

(c4 66) subst ( ww, p4al A 2+p3a2al A 2 , %) ; 

2 2 2 
p3a2al 2 cos (2 gh) p3a2al p4al sin (2 gh) p4al (p3a2al - p4al ) 

(d4 66) + + 

ww 2 2 

ww ww 

(c467) fn: subst (wwm, p3a2al A 2-p4al A 2 , %) ; 

2 

sin (2 gh) p4al wwm p3a2al 2 cos (2 gh) p3a2al p4al 

(d4 67 ) + + 

2 ww 2 

ww ww 

(c468) fn : combine (fn) ; 

2 

sin (2 gh) p4al wwm + 2 cos (2 gh) p3a2al p4al p3a2al 


(d4 68 ) + 

2 ww 

ww 


(c469) /* 4 * k(12,22) 2 to 1 term f22 gsyy *** c21f22y ************************* 
kl222 : (a2 / (4*%pi*all) ) *fn; 

2 

sin (2 gh) p4al wwm + 2 cos (2 gh) p3a2al p4al p3a2al 

a 2 ( + } 

2 ww 

ww 

(d469) 

4 %pi all 

(c470) /* where 

p3a2al=-p3+a2*gt-al*gc*cos (gh) 

p4al=p4-al*gc*sin (gh) 

ww=p4al A 2+p3a2al A 2 

wwm=p3a2al A 2-p4al A 2 
*/ 

fortran (p3a2al : -p3+a2*gt ( j) -al*gc (i) *dcos (gh) ) $ 
p3a2al = -p3+a2*gt ( j) -al*dcos (gh) *gc (i) 

(c471) fortran (p4al :p4-al*gc (i) *dsin (gh) ) $ 
p4al = p4-al*dsin (gh) *gc (i) 

(c472 ) fortran (ww :'p4al^2+'p3a2al A 2)$ 
ww = p4al**2+p3a2al**2 


(c473) fortran (wwm: 'p3a2al A 2-'p4al A 2) $ 
wwm = p3a2al**2-p4al**2 



(c4 74) subst (dcos, cos, first (first (fn) )) $ 

(c475) subst (dsin, sin, %) $ 

(c476) fortran (parti :%) $ 

parti = dsin (2*gh) *p4al*wwm+2*dcos (2*gh) *p3a2al*p4al**2 

(c4 77 ) fortran (parti : 'parti / ( ' ww / '2) ) $ 
parti = partl/ww**2 

(c478) subst (dcos, cos, last (fn) ) $ 

(c479) subst (dsin, sin, %) $ 

(c480) fortran (part2 :%) $ 
part2 = p3a2al/ww 

(c481 ) kill (xl ) $ 

(c482 ) kill (x2 ) $ 

(c483 ) kill (yl ) $ 

(c484) kill (y2) $ 

(c485) kill (t 1 ) $ 

(c486) kill (t2) $ 

(c487) kill (p3a2al) $ 

(c488) kill (p4al) $ 

(c489) kill (ww) $ 

(c490) kill (wwm) $ 

(c491) x2:a2*gc; 

(d491) a2 gc 

(c492) tl:al*gt; 

(d492) al gt 

(c493) t2:a2*gt; 

(d493) a2 gt 

(c494) gsxl : ' ' gsxxpfltl; 

2 2 

2 yl yl ( (al gt - xl) - yl ) 

(d4 94 ) + 

2 2 2 2 2 
yl + (al gt - xl) (yl + (al gt - xl) ) 

(c495) gsyl : ' ' gsyypf ltl; 

2 2 

yl ( (al gt - xl) - yl ) 


2 2 2 
(yl + (al gt - xl) ) 


(d495) 



(c496) gtxyl : ' ' gsxypf ltl; 


(d4 96 ) 


2 (al gt - xl) yl 


al gt - xl 


2 2 2 2 2 
(yl + (al gt - xl) ) yl + (al gt - xl) 


(c497) xl :pl+x2*cos (gh) -y2*sin (gh) ; 

(d497) - sin(gh) y2 + pi + a2 gc cos(gh) 

(c498) yl :p2+x2*sin (gh) +y2*cos (gh) ; 

(d498) cos (gh) y2 + p2 + a2 gc sin (gh) 

(c499) /* 5 * k(21,ll) 1 to 2 term fll gtxy *** cl2fllxy *************************** 
gtxyl 2b : - ( ' gsxl-' gsyl ) /2*sin (2*gh) +' gtxyl *cos (2*gh) ; 

sin(2 gh) (gsxl - gsyl) 

2 

(c500) gtxyl2 : - ( ' ' gsxl-' ' gsyl) /2*sin (2*gh) + ' ' gtxyl* cos (2*gh) ; 

2 

(d500) cos (2 gh) (2 (cos(gh) y2 + p2 + a2 gc sin(gh)) 

(sin(gh) y2 - pi + al gt - a2 gc cos (gh) ) 


(d499) 


cos (2 gh) gtxyl - 


2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) , 2) 

- (sin(gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

/((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) )) 

- sin(2 gh) (2 (cos(gh) y2 + p2 + a2 gc sin(gh)) 

2 

/((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) ) + 2 (cos (gh) y2 + p2 + a2 gc sin (gh) ) 

2 

((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

- (cos (gh) y2 + p2 + a2 gc sin(gh)) ) 


2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin (gh) ) , 2 ) ) / 2 
(c501) y2 : 0 ; 



(d501) 


0 


(c502) ev(gtxyl2); 


(d502 ) cos (2 gh) (- 


2 (- pi + al gt - a2 gc cos (gh) ) (p2 + a2 gc sin(gh)) 

2 2 2 
( (p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) ) 


- pi + al gt - a2 gc cos(gh) 


2 2 
(p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) 


- sin(2 gh) (- 


2 (p2 + a2 gc sin(gh)) 


2 2 
(p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos (gh) ) 


+ 2 (p2 + a2 gc sin(gh)) ((-pi + al gt - a2 gc cos(gh) 


- (p2 + a2 gc sin(gh)) ) 

2 2 2 
/ ( (p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos (gh) ) ) )/2 

(c503) subst (plala2, -pl+al*gt-a2*gc*cos (gh) , %) ; 


(d503) cos (2 gh) (- 


2 plala2 (p2 + a2 gc sin(gh)) 

2 2 2 
( (p2 + a2 gc sin(gh)) + plala2 ) 


plala2 


(- 


2 2 

(p2 + a2 gc sin(gh)) + plala2 
2 (p2 + a2 gc sin(gh)) 

2 2 

(p2 + a2 gc sin(gh)) + plala2 


■) - sin (2 gh) 


2 2 
2 (p2 + a2 gc sin(gh)) (plala2 - (p2 + a2 gc sin(gh)) ) 

+ )/2 

2 2 2 

( (p2 + a2 gc sin(gh)) + plala2 ) 

(c504) subst (p2a2,p2+a2*gc*sin (gh) , %) ; 

2 


(d504) cos(2 gh) (• 


2 plala2 p2a2 


plala2 


-) 


2 2 2 2 2 
(p2a2 + plala2 ) p2a2 + plala2 


sin(2 gh) (- 


2 p2a2 2 p2a2 (plala2 - p2a2 ) 

+ } 

2 2 2 2 2 



p2a2 + plala2 


(p2a2 + plala2 ) 


(c505) 
(d505 ) 

(c506) 

(d506) 
(c507 ) 
(d507 ) 
(c508 ) 
(d508 ) 
(c509) 

(d509) 
(c510 ) 
(d510) 

(c511) 


2 


subst ( ww, p2a2 A 2+plala2 A 2 , % ) ; 

2 


2 plala2 p2a2 plala2 

cos (2 gh) ( ) 

2 ww 

ww 


sin (2 gh) 


2 p2a2 2 p2a2 (plala2 


( + 

ww 2 

ww 


2 

p2a2 ) 
) 


f n : subst (wwm, plala2 A 2-p2a2 A 2 , % ) ; 


cos (2 gh) 


2 

2 plala2 p2a2 

( 

2 

ww 


plala2 

) 

ww 


2 p2a2 wwm 2 p2a2 

sin(2 gh) ( + ) 

2 ww 

ww 


2 


first (f n) ; 

2 


2 plala2 p2a2 plala2 

cos (2 gh) ( ) 

2 ww 

ww 


n2 : mult thru (%) ; 

2 

2 cos (2 gh) plala2 p2a2 cos (2 gh) plala2 


2 ww 

ww 


rest (f n, 1) ; 

2 p2a2 wwm 2 p2a2 

sin (2 gh) < + ) 

2 ww 

ww 

2 

n3 :multthru {%) ; 

sin (2 gh) p2a2 wwm sin (2 gh) p2a2 

2 ww 

ww 


n : n2+n3; 


(d511) 


sin (2 gh) p2a2 wwm sin (2 gh) p2a2 cos (2 gh) plala2 



2 


ww 


ww 


ww 


2 

2 cos (2 gh) plala2 p2a2 



2 

ww 


(c5l2) 


(d5l2) 


f n : combine (n) ; 

2 

2 cos (2 gh) plala2 p2a2 - sin (2 gh) p2a2 wwm 


2 

ww 


- sin (2 gh) p2a2 - cos (2 gh) plala2 


(c5l3) /* 5 * k(21,ll) 1 to 2 term fll gtxy *** cl2fllxy *************************** 
k2111 : (al/ (4*%pi*all) ) *fn; 


2 

2 cos (2 gh) plala2 p2a2 - sin (2 gh) p2a2 wwm 

(d513) al ( 

2 

ww 

- sin (2 gh) p2a2 - cos (2 gh) plala2 

+ ) / (4 %pi all) 

ww 

(c514) /* where 

plala2=-pl+al*gt-a2*gc*cos (gh) 

p2a2=p2+a2*gc*sin (gh) 

ww=p2a2 A 2+plala2 A 2 

wwm=plala2 A 2-p2a2 A 2 
*/ 

fortran (plala2 : -pl+al*gt ( j ) -a2*gc (i) *dcos (gh) ) $ 
plala2 = -pl+al*gt ( j) -a2*dcos (gh) *gc (i) 

(c515) fortran (p2a2 :p2+a2*gc (i) *dsin (gh) ) $ 
p2a2 = p2+a2*dsin (gh) *gc ( i) 

(c516) fortran (ww : ' p2a2 A 2+' plala2 A 2 ) $ 
ww = p2a2**2+plala2**2 

(c517) fortran (wwm: f plala2 A 2-'p2a2 A 2) $ 
wwm = plala2**2-p2a2**2 

(c518) subst (dcos, cos, first (first (fn) ) ) $ 

(c519) subst (dsin, sin, %) $ 

(c520) fortran (parti :%) $ 

parti = 2*dcos (2*gh) *plala2*p2a2**2-dsin(2*gh) *p2a2*wwm 


(c521) fortran (parti : 'parti/ ( ' ww A 2) ) $ 
parti = partl/ww**2 



(c522) subst (dcos, cos, last (fn) ) $ 
(c523) subst (dsin, sin, %) $ 


(c524 ) 

fortran (part2 : %) $ 

part2 = (-dsin (2*gh) *p2a2-dcos (2*gh) *plala2 ) / ww 

(c525) 

kill (xl) $ 



( c526 ) 

kill (x2 ) $ 



(c527 ) 

kill (yl) $ 



(c528) 

kill (y2) $ 



(c529) 

kill (tl) $ 



(c530 ) 

kill (t 2 ) $ 



(c531) 

kill (plala2) $ 



( c532 ) 

kill (p2a2) $ 



(c533) 

kill (ww) $ 



(c534) 

kill (wwm) $ 



(c535) 

x2 : a2*gc; 



(d535) 



a2 gc 

(c536) 

tl :al*gt; 



(d536) 



al gt 

(c537 ) 

t2 :a2*gt; 



(d537 ) 



a2 gt 

(c538 ) 

gsxl : ' ' gsxxpf 2t 1 ; 



(d538 ) 

al gt - xl 


2 (al gt - xl) 

2 

yl + (al gt - xl) 

2 

2 

(yl + (al gt - 

(c539 ) 

gsyl : ' ' gsyypf2tl; 



(d539) 

al gt - xl 


2 (al gt - xl) 

2 

yl + (al gt - xl) 

2 

2 

(yl + (al gt - 

(c540 ) 

gtxyl : ' ' gsxypf 2tl; 



(d540 ) 

yl ( (al 

gt 

2 2 
- xl) - yl ) 

2 

(yl + 

2 2 

(al gt - xl) ) 

(c541) 

xl :pl+x2*cos (gh) -y2*sin (gh) ; 




2 

yi 


2 2 
xl) ) 


2 

yi 


2 2 
xl) ) 



(d541) 


- sin(gh) y2 + pi + a2 gc cos(gh) 


(c542 ) yl :p2+x2*sin (gh) +y2*cos (gh) ; 

(d542) cos (gh) y2 + p2 + a2 gc sin (gh) 

(c543) /* 6 * k(21,12) 1 to 2 term f!2 gtxy *** cl2f!2xy ************** 

gtxyl2b : - ( ' gsxl-' gsyl ) /2*sin (2*gh) + ' gtxyl *cos (2*gh) ; 

sin(2 gh) (gsxl - gsyl) 

(d543) cos (2 gh) gtxyl 

2 

(c544) gtxyl2 : - { ' ' gsxl-' ' gsyl ) /2*sin ( 2*gh) + ' ' gtxyl* cos (2*gh) ; 

(d544) cos (2 gh) (cos (gh) y2 + p2 + a2 gc sin(gh)) 

2 

((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

- (cos (gh) y2 + p2 + a2 gc sin(gh)) ) 

2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) , 2) 

2 

+ 2 sin (2 gh) (cos (gh) y2 + p2 + a2 gc sin(gh)) 

(sin(gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) , 2) 

(c545) y2:0; 

(d545) 0 

(c546) ev(gtxyl2); 

2 

(d546) cos (2 gh) (p2 + a2 gc sin(gh)) ((-pi + al gt - a2 gc cos(gh)) 

2 

- (p2 + a2 gc sin(gh)) ) 

2 2 2 
/ ( (p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos (gh) ) ) 

2 

2 sin (2 gh) (- pi + al gt - a2 gc cos(gh)) (p2 + a2 gc sin(gh)) 

+ 

2 2 2 
( (p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) ) 

(c547 ) subst (plala2 , -pl+al*gt-a2*gc*cos (gh) , % ) ; 

2 2 
cos (2 gh) (p2 + a2 gc sin(gh)) (plala2 - (p2 + a2 gc sin(gh)) ) 



(d547 ) 


2 2 2 

( (p2 + a2 gc sin(gh)) + plala2 ) 

2 

2 sin (2 gh) plala2 (p2 + a2 gc sin(gh)) 


2 2 2 

( (p2 + a2 gc sin(gh)) + plala2 ) 


(c548) 
(d548) 
(c549) 
(d549) 
(c550 ) 
(d550 ) 


subst (p2a2 , p2+a2*gc*sin (gh) , %) ; 

2 2 2 

2 sin (2 gh) plala2 p2a2 cos (2 gh) p2a2 (plala2 - p2a2 ) 

+ 

2 2 2 2 2 2 
(p2a2 + plala2 ) (p2a2 + plala2 ) 

subst (ww,p2a2 A 2+plala2 A 2, %) ; 

2 2 2 

2 sin (2 gh) plala2 p2a2 cos (2 gh) p2a2 (plala2 - p2a2 ) 

+ 

2 2 
ww ww 

f n : subst (wwm, plala2 A 2-p2a2 A 2 , %) ; 

2 

cos (2 gh) p2a2 wwm 2 sin (2 gh) plala2 p2a2 

+ 

2 2 
ww ww 


(c55l) fn : combine (fn) ; 


(d551 ) 


2 

cos (2 gh) p2a2 wwm + 2 sin (2 gh) plala2 p2a2 


2 

ww 


(c552) /* 6 * k(21,12) 1 to 2 term fl2 gtxy *** cl2fl2xy *********************** 
k2112 : (al/ (4*%pi*all) ) *fn; 


(d552 ) 


2 

al (cos (2 gh) p2a2 wwm + 2 sin (2 gh) plala2 p2a2 ) 


2 

4 %pi all ww 


(c553) /* where 

plala2=-pl+al*gt-a2*gc*cos (gh) 

p2a2=p2+a2*gc*sin (gh) 

ww=p2a2 A 2+plala2 A 2 

wwm=plala2 A 2-p2a2 A 2 
*/ 


fortran (plala2 : -pl+al*gt ( j) -a2*gc (i) *dcos (gh) ) $ 
plala2 = -pl+al*gt ( j) -a2*dcos (gh) *gc (i) 

(c554 ) fortran (p2a2 :p2+a2*gc (i) *dsin (gh) ) $ 



p2a2 = p2+a2*dsin (gh) *gc (i) 

(c555 ) fortran (ww : ' p2a2 A 2+' plala2 A 2 ) $ 
ww = p2a2**2+plala2**2 

(c556 ) fortran ( wwm: 'plala2 A 2-'p2a2 A 2) $ 
wwm = plala2**2-p2a2**2 

(c557 ) subst (dcos, cos, first (fn) ) $ 

(c558) subst (dsin, sin, %) $ 

(c559) fortran (parti :%) $ 

parti = dcos (2*gh) *p2a2*wwm+2*dsin (2*gh) *plala2*p2a2**2 

(c560) fortran (parti : 'parti/ ( ' ww A 2) ) $ 
parti = partl/ww**2 

(c561) kill (xl ) $ 

(c562) kill (x2) $ 

(c563) kill (yl) $ 

(c564 ) kill (y2) $ 

(c565) kill (tl) $ 

(c566 ) kill (t2) $ 

(c567 ) kill (plala2) $ 

(c568 ) kill (p2a2 ) $ 

(c569 ) kill (ww) $ 

(c570) kill (wwm) $ 

(c571) x2:a2*gc; 

(d571) a2 gc 

(c572) tl:al*gt; 

(d572 ) al gt 

(c573) t2:a2*gt; 

(d573) a2 gt 

(c574) gsxl : ' ' gsxxpf ltl; 


(d574 ) 


2 2 

2 yl yl ( (al gt - xl) - yl ) 


2 2 2 2 2 
yl + (al gt - xl) (yl + (al gt - xl) ) 


(c575) gsyl : ' ' gsyypf ltl; 

2 2 

yl ( (al gt - xl) - yl ) 


2 2 2 
(yl + (al gt - xl) ) 


(d575) 



(c576) gtxyl : ' ' gsxypf ltl; 


2 


2 (al gt - xl) yl al gt - xl 

( d 576 ) 


2 2 2 2 2 
(yl + (al gt - xl) ) yl + (al gt - xl) 

(c577) xl :pl+x2*cos (gh) -y2*sin (gh) ; 

(d577) - sin(gh) y2 + pi + a2 gc cos(gh) 

(c578 ) yl :p2+x2*sin (gh) +y2*cos (gh) ; 

(d578) cos (gh) y2 + p2 + a2 gc sin (gh) 

(c579) /* 7 * k(22,ll) 1 to 2 term fll gsyy *** cl2flly *************************** 

gsyl2b : ( ' gsxl+' gsyl ) /2- ( ' gsxl-' gsyl) /2*cos (2*gh) gtxyl *sin (2*gh) / 

gsyl + gsxl cos (2 gh) (gsxl - gsyl) 

(d579) - sin(2 gh) gtxyl + 

2 2 

(c580) gsyl2: ( ' ' gsxl+' ' gsyl) /2- (' ' gsxl-' 'gsyl ) /2*cos (2*gh) gtxyl *sin (2*gh) ; 


2 


(d580) - sin (2 gh) (2 (cos(gh) y2 + p2 + a2 gc sin(gh)) 
(sin(gh) y2 - pi + al gt - a2 gc cos (gh) ) 


2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) , 2) 

- (sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

/((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 


2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) )) 

- cos (2 gh) (2 (cos(gh) y2 + p2 + a2 gc sin(gh)) 


2 

/((sin(gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) ) + 2 (cos (gh) y2 + p2 + a2 gc sin(gh)) 

2 

((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

- (cos (gh) y2 + p2 + a2 gc sin(gh)) ) 


2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin (gh) ) , 2 ) ) / 2 
+ (cos (gh) y2 + p2 + a2 gc sin (gh) ) 



2 

/((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) ) 

(c581) y2:0; 

(d581 ) 0 

(c582) ev(gsyl2); 

2 (p2 + a2 gc sin (gh) ) 

(d582) - cos(2 gh) ( 

2 2 
(p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) 

2 

+ 2 (p2 + a2 gc sin(gh)) ((-pi + al gt - a2 gc cos(gh)) 

2 

- (p2 + a2 gc sin(gh)) ) 

2 2 2 
/ ( (p2 + a2 gc sin (gh) ) + (- pi + al gt - a2 gc cos (gh) ) ) )/2 

2 

2 (- pi + al gt - a2 gc cos (gh) ) (p2 + a2 gc sin(gh)) 

- sin(2 gh) ( 

2 2 2 
( (p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) ) 

- pi + al gt - a2 gc cos(gh) 

) 

2 2 

(p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) 

p2 + a2 gc sin(gh) 

+ 

2 2 

(p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) 

(c583) subst (plala2, -pl+al*gt-a2*gc*cos (gh) , %) ; 

2 (p2 + a2 gc sin(gh)) 

(d583) - cos(2 gh) ( 

2 2 

(p2 + a2 gc sin (gh) ) + plala2 

2 2 

2 (p2 + a2 gc sin(gh)) (plala2 - (p 2 + a2 gc sin(gh)) ) 

+ )/2 

2 2 2 

( (p2 + a2 gc sin(gh)) + plala2 ) 

2 

2 plala2 (p2 + a2 gc sin(gh)) 

- sin(2 gh) ( 

2 2 2 

( (p2 + a2 gc sin(gh)) + plala2 ) 

plala2 p2 + a2 gc sin(gh) 

) + 

2 2 2 2 
(p2 + a2 gc sin(gh)) + plala2 (p2 + a2 gc sin(gh)) + plala2 



(c584 ) subst (p2a2 ,p2+a2*gc*sin (gh) , %) ; 


(d584 ) 

(c585) 

(d585) 

(c586) 

(d586) 

( c587 ) 
(d587) 
(c588 ) 

(d588 ) 
(c589 ) 


cos (2 gh) 


2 2 

2 p2a2 2 p2a2 (plala2 - p2a2 ) 

( + ) 

2 2 2 2 2 
p2a2 + plala2 (p2a2 + plala2 ) 


2 


2 

2 plala2 p2a2 plala2 

- sin (2 gh) ( ) 

2 2 2 2 2 


(p2a2 + plala2 ) p2a2 + plala2 

subst (ww, p2a2 / '2+plala2 A 2 , %) ; 


cos (2 gh) 


2 2 

2 p2a2 2 p2a2 (plala2 - p2a2 ) 

( + ) 

ww 2 


ww 


p2a2 

2 2 
p2a2 + plala2 


2 

2 


2 plala2 p2a2 plala2 p2a2 

- sin(2 gh) { ) + 

2 ww ww 

ww 


fn : subst (wwm, plala2 A 2-p2a2 A 2 , %) ; 

2 p2a2 wwm 2 p2a2 

cos (2 gh) ( + ) 

2 ww 2 


ww 2 plala2 p2a2 plala2 

sin (2 gh) ( ) 

2 2 ww 

ww 


p2a2 

ww 

nl .-last (fn) ; 

p2a2 

ww 

first (fn) ; 

2 p2a2 wwm 2 p2a2 

cos (2 gh) ( + ) 

2 ww 

ww 

2 

n2 :multthru (%) ; 


(d589 ) 


cos (2 gh) p2a2 wwm cos (2 gh) p2a2 



2 


ww 


ww 

(c590) rest (f n, 1) ; 

2 

p2a2 2 plala2 p2a2 plala2 


( d5 90) sin (2 gh) ( ) 

ww 2 ww 


ww 

(c591) last (%) ; 

2 


2 plala2 p2a2 plala2 

(d591) - sin (2 gh) ( ) 

2 ww 

ww 


(c592) n3 :multthru (%) ; 

2 

sin (2 gh) plala2 2 sin (2 gh) plala2 p2a2 

(d592) 

ww 2 

ww 

(c593) n:nl+n2+n3; 

cos (2 gh) p2a2 wwm cos (2 gh) p2a2 p2a2 sin (2 gh) plala2 


(d593) + + 

2 ww ww ww 

ww 


2 

2 sin (2 gh) plala2 p2a2 


ww 

(c594) n : combine (n) ; 

2 

- cos (2 gh) p2a2 wwm - 2 sin (2 gh) plala2 p2a2 

(d594) 

2 

ww 

- cos (2 gh) p2a2 + p2a2 + sin (2 gh) plala2 


(c595) tl:first(%); 

2 

- cos (2 gh) p2a2 wwm - 2 sin (2 gh) plala2 p2a2 

(d595) 

2 

ww 

(c596) last (n) ; 

- cos (2 gh) p2a2 + p2a2 + sin (2 gh) plala2 

( d596) 

ww 



(c597) t2 :part (%, 1) ; 


(d597) - cos (2 gh) p2a2 + p2a2 + sin (2 gh) plala2 

(c598) first (%) +f irst (rest (%, 1) ) ; 

(d598) p2a2 - cos (2 gh) p2a2 

(c599) factor (%) ; 

(d599 ) - (cos (2 gh) - 1) p2a2 

(c600) trigexpand (%) ; 

2 2 

(d600) - (- sin (gh) + cos (gh) - 1) p2a2 

(c601 ) t2 : (trigsimp (% ) +last (t2) )/ww; 

2 

2 sin (gh) p2a2 + sin (2 gh) plala2 

(d601) 

ww 

(c602) /* 7 * k(22,ll) 1 to 2 term fll gsyy *** c!2flly **************************** 
k2211 : (al/(4*%pi*all) )*(tl+t2); 


(d602 ) al ( 


2 

- cos (2 gh) p2a2 wwm - 2 sin (2 gh) plala2 p2a2 


2 

ww 


2 

2 sin (gh) p2a2 + sin (2 gh) plala2 

+ ) / (4 %pi all) 

ww 

(c603) /* where 

plala2=-pl+al*gt-a2*gc*cos (gh) 

p2a2=p2+a2*gc*si n (gh) 

ww=p2a2 A 2+plala2 A 2 

wwm=plala2 A 2-p2a2 A 2 
*/ 

fortran (plala2 : -pl+al*gt ( j ) -a2*gc (i) *dcos (gh) ) $ 
plala2 = -pl+al*gt ( j ) -a2*dcos (gh) *gc (i) 

(c604) fortran (p2a2 :p2+a2*gc (i) *dsin (gh) ) $ 
p2a2 = p2+a2*dsin (gh) *gc (i) 

(c605) fortran (ww : ' p2a2 A 2 + f plala2 A 2 ) $ 
ww = p2a2**2+plala2**2 

(c606) fortran (wwm: , plala2 A 2- , p2a2 A 2) $ 
wwm = plala2**2-p2a2**2 

(c607) subst (dcos, cos, first (tl) ) $ 

(c608) subst (dsin, sin, %) $ 

(c609) fortran (parti :%) $ 



parti = -dcos (2*gh) *p2a2*wwm-2*dsin (2*gh) *plala2*p2a2**2 


(c610) fortran (parti : 'parti/ ( ' ww A 2) ) $ 
parti = partl/ww**2 

(c611) subst (dcos, cos, t2) $ 

(c612) subst (dsin, sin, %) $ 

(c613) fortran (part2 :%) $ 

part2 = (2*dsin (gh) **2*p2a2+dsin (2*gh) *plala2 ) /ww 


(c614 ) 

kill (xl) $ 






(c615) 

kill (x2 ) $ 






(c616) 

kill (yl) $ 






(c617) 

kill (y2 ) $ 






(c618) 

kill (tl) $ 






( c 6 19 ) 

kill (t2) $ 






(c620 ) 

kill (plala2 ) $ 






(c621) 

kill (p2a2) $ 






(c622 ) 

kill (ww) $ 






(c623 ) 

kill (wwm) $ 






(c624 ) 

x2 :a2*gc; 






( d6 24) 





a2 gc 


(c625) 

tl :al*gt; 






(d625 ) 





al gt 


(c626) 

t2 :a2*gt; 






(d626) 





a2 gt 


(c627 ) 

gsxl : ' ' gsxxpf2tl; 





(d627 ) 


al gt 

- xl 


2 (al 

2 

gt - xl) yl 

2 

yi 

+ (al 

gt - 

2 

xl) 

2 

(yl + 

2 2 

(al gt - xl) ) 

(c628 ) 

gsyl : ' ' gsyypf2tl; 





(d628 ) 


al gt 

- xl 


2 (al 

2 

gt - xl) yl 

2 

yi 

+ (al 

gt - 

2 

xl) 

2 

(yl + 

2 2 

(al gt - xl) ) 

(c629) 

gtxyl : ' ' gsxypf 2tl; 





(d629) 



yl ( (al gt 

2 

- xl) - 

2 

yi ) 


2 2 



(yl + (al gt - xl) ) 


(c6 30) xl :pl+x2*cos (gh) -y2*sin (gh) ; 

(d630) - sin(gh) y2 + pi + a2 gc cos(gh) 

(c63l) yl :p2+x2*sin (gh) +y2*cos (gh) ; 

(d63l) cos (gh) y2 + p2 + a2 gc sin (gh) 

(c632) /* 8 * k(22,12) 1 to 2 term fl2 gsyy *** cl2fl2y **************************** 

gsyl2b : { ' gsxl+' gsyl) /2- ( ' gsxl-' gsyl) /2*cos (2*gh) -' gtxyl*sin (2*gh) ; 

gsyl + gsxl cos (2 gh) (gsxl - gsyl) 

(d632) - sin (2 gh) gtxyl + 

2 2 


(c633) gsyl 2 : ( ' ' gsxl+' ' gsyl) /2- ( ' 'gsxl-' ' gsyl) /2*cos (2*gh) -' ' gtxyl* sin (2*gh) ; 
(d633) (sin(gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

/((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) ) - sin (2 gh) 


2 

(cos (gh) y2 + p2 + a2 gc sin(gh)) ((sin(gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

- (cos (gh) y2 + p2 + a2 gc sin(gh)) ) 

2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos(gh)) 

2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) , 2) 


2 

+ 2 cos (2 gh) (cos (gh) y2 + p2 + a2 gc sin(gh)) 
(sin(gh) y2 - pi + al gt - a2 gc cos (gh) ) 

2 

/expt ( (sin (gh) y2 - pi + al gt - a2 gc cos (gh) ) 


2 

+ (cos (gh) y2 + p2 + a2 gc sin(gh)) , 2) 
(c634 ) y2 : 0 ; 

(d634) 0 

(c635) ev (gsyl2) ; 


(d635) 


- pi + al gt - a2 gc cos(gh) 


2 2 
(p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) 


2 

- sin (2 gh) (p2 + a2 gc sin(gh)) ((-pi + al gt - a2 gc cos(gh)) 


2 



(p2 + a2 gc sin(gh)) ) 


2 2 2 
/ ( (p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos (gh )) ) 

2 

2 cos {2 gh) (- pi + al gt - a2 gc cos(gh)) (p2 + a2 gc sin (gh) ) 

+ 

2 2 2 
( (p2 + a2 gc sin(gh)) + (- pi + al gt - a2 gc cos(gh)) ) 

(c636) subst (plala2, -pl+al*gt-a2*gc*cos (gh) , %) ; 

plala2 

(d636 ) 

2 2 

(p2 + a2 gc sin(gh)) + plala2 

2 2 
sin (2 gh) (p2 + a2 gc sin(gh)) (plala2 - (p2 + a2 gc sin(gh)) ) 

2 2 2 

( (p2 + a2 gc sin(gh)) + plala2 ) 

2 

2 cos (2 gh) plala2 (p2 + a2 gc sin(gh)) 

+ 

2 2 2 

( (p2 + a2 gc sin(gh)) + plala2 ) 

(c637) subst (p2a2,p2+a2*gc*sin (gh) , %) ; 

2 

plala2 2 cos (2 gh) plala2 p2a2 

(cl637) + 

2 2 2 2 2 
p2a2 + plala2 (p2a2 + plala2 ) 


2 2 

sin (2 gh) p2a2 (plala2 - p2a2 ) 


2 2 2 
(p2a2 + plala2 ) 

(c638 ) subst (ww, p2a2 A 2+plala2 A 2 , % ) ; 

2 2 2 
plala2 2 cos (2 gh) plala2 p2a2 sin (2 gh) p2a2 (plala2 - p2a2 ) 


(d638) + 

ww 2 2 

ww ww 


(c639 ) f n : subst (wwm, plala2 A 2-p2a2 A 2 , % ) ; 

2 

sin (2 gh) p2a2 wwm plala2 2 cos (2 gh) plala2 p2a2 

(d639) + + 

2 ww 2 

ww ww 

(c640) f n : combine (fn) ; 

2 

2 cos (2 gh) plala2 p2a2 - sin (2 gh) p2a2 wwm plala2 


(d640) + 

2 ww 



ww 


(c641) /* 8 * k(22,12) 1 to 2 term fl2 gsyy *** cl2f!2y **************************** 
k2212 : (al/ (4*%pi*all) ) *fn; 


(d641 ) 


al 


2 

2 cos (2 gh) plala2 p2a2 


ww 


sin (2 gh) p2a2 wwm plala2 

+ ) 

ww 


4 %pi all 


(c642) /* where 

plala2=-pl+al*gt-a2*gc*cos (gh) 


p2a2=p2+a2*gc*sin (gh) 


ww=p2a2 A 2+plala2 A 2 

wwm=plala2 A 2-p2a2 A 2 
*/ 


fortran (plala2 : -pl+al*gt ( j) -a2*gc (i) *dcos (gh) ) $ 
plala2 = -pl+al*gt ( j ) -a2*dcos (gh) *gc (i) 

(c643) fortran (p2a2 :p2+a2*gc (i) *dsin (gh) ) $ 
p2a2 = p2+a2*dsin (gh) *gc (i) 

(c64 4) fortran (ww: ' p2a2 A 2 + ' plala2 A 2 ) $ 
ww = p2a2**2+plala2**2 

(c645) fortran (wwm: , plala2 A 2- , p2a2 A 2) $ 
wwm = plala2**2-p2a2**2 

(c64 6) subst (dcos, cos, first (f irst (fn) ) ) $ 

(c647) subst (dsin, sin, %) $ 

(c648) fortran (parti :%) $ 

parti = 2*dcos (2*gh) *plala2*p2a2**2-dsin (2*gh) *p2a2*wwm 

(c64 9) fortran (parti : 'parti/ ( ' ww A 2) ) $ 
parti = partl/ww**2 

(c650) subst (dcos, cos, last (fn) ) $ 

(c651) subst (dsin, sin, %) $ 

(c652) fortran (part2 :%) $ 
part2 = plala2/ww 

(c653) kill (xl) $ 

(c654 ) kill (x2) $ 

(c655) kill (yl) $ 

(c656) kill (y2) $ 

(c657) kill (tl) $ 

(c658 ) kill (t2) $ 

(c659 ) kill (plala2) $ 



(c660) kill (p2a2) $ 

(c661 > kill (ww) $ 

(c 662) kill (wwm) $ 

(c663) closef ile ( totallp) 



